Introduction

My name is Dustin Bracy, and this is a time series analysis of COVID-19 data in Louisiana and the United States spanning an 8 month period from the declaration of the pandemic on March 11, 2020 through November 11, 2020. Below I will explore the difference between new daily positive cases and positivity rates, and why one may be more advantageous. Next I’ll walk through univiariate and multivariate analyses of both Louisiana and national numbers. Finally I’ll wrap up with a comparison of model performance for each use case.

Acquire and Load the Data

The COVID-19 data used here comes from two primary sources: the Louisiana Department of Health Coronavirus website, and The Covid Tracking Project website. Additional temperature, windspeed and other related data for Louisiana was collected manually from historical data for the MSY (New Orleans regional airport) section of Weather Underground.

The data comes in pretty good shape, with only a handful of missing values in the US dataset. Fortunately, there are no missing values in response variables of interest (i.e. the total number of daily tests and the number of daily positive cases).

The data ingestion pipeline is fairly simple: download the data directly from the api or website, remove data outside of the March 11 to November 11 window, and add positivity metrics.

# Get latest US Data & sort by date ascending
download.file('https://covidtracking.com/data/download/national-history.csv',
              destfile = './data/national-history.csv')

us <- read_csv('./data/national-history.csv')
## 
## ── Column specification ─────────────────────────────────────────────────────
## cols(
##   date = col_date(format = ""),
##   death = col_double(),
##   deathIncrease = col_double(),
##   inIcuCumulative = col_double(),
##   inIcuCurrently = col_double(),
##   hospitalizedIncrease = col_double(),
##   hospitalizedCurrently = col_double(),
##   hospitalizedCumulative = col_double(),
##   negative = col_double(),
##   negativeIncrease = col_double(),
##   onVentilatorCumulative = col_double(),
##   onVentilatorCurrently = col_double(),
##   positive = col_double(),
##   positiveIncrease = col_double(),
##   recovered = col_double(),
##   states = col_double(),
##   totalTestResults = col_double(),
##   totalTestResultsIncrease = col_double()
## )
us <- us[order(us$date),]


# Get latest LA Data & consolidate to state level
download.file('https://ldh.la.gov/assets/oph/Coronavirus/data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx',
              destfile = './data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx')

la <- readxl::read_xlsx('./data/LA_COVID_TESTBYDAY_PARISH_PUBLICUSE.xlsx')
la <- la %>%
  group_by(date = `Lab Collection Date`) %>%
  summarise(
    tests = sum(`Daily Test Count`),
    positive = sum(`Daily Positive Test Count`),
    negative = sum(`Daily Negative Test Count`),
    cases = sum(`Daily Case Count`)
  )
## `summarise()` ungrouping output (override with `.groups` argument)
# Create positivity features
la$positivity = la$positive / la$tests
us$positivity <- us$positiveIncrease / us$totalTestResultsIncrease


# Get temp data from WeatherUnderground.com (historical data for MSY)
## Note Nov8 data missing from site, so the imputed average of Nov7+9 was taken
xdf <- read_csv('./data/temp.csv')
## 
## ── Column specification ─────────────────────────────────────────────────────
## cols(
##   Date = col_character(),
##   MaxTempF = col_double(),
##   AvgTempF = col_double(),
##   MinTempF = col_double(),
##   MaxDewPointF = col_double(),
##   AvgDewPointF = col_double(),
##   MinDewPointF = col_double(),
##   MaxHumidityPercent = col_double(),
##   AvgHumidityPercent = col_double(),
##   MinHumidityPercent = col_double(),
##   MaxWindSpeedMPH = col_double(),
##   AvgWindSpeedMPH = col_double(),
##   MinWindSpeedMPH = col_double(),
##   MaxPressureHg = col_double(),
##   AvgPressureHg = col_double(),
##   MinPressureHg = col_double(),
##   PrecipitationIN = col_double(),
##   Weekday = col_character(),
##   Gathering = col_character()
## )
xdf$GatheringBinary <- as.numeric(ifelse(is.na(xdf$Gathering), 0,1))

# Trim all datasets to span March 11, 2020 through Nov 11, 2020
la <- la[11:256,]
us <- us[50:295,]
xdf <- xdf[11:256,]

# Check for missing values
la %>% is.na() %>% summary()
##     date           tests          positive        negative      
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246       FALSE:246       FALSE:246      
##    cases         positivity     
##  Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246
us %>% is.na() %>% summary()
##     date           death         deathIncrease   inIcuCumulative
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246       FALSE:246       FALSE:232      
##                                                  TRUE :14       
##  inIcuCurrently  hospitalizedIncrease hospitalizedCurrently
##  Mode :logical   Mode :logical        Mode :logical        
##  FALSE:231       FALSE:246            FALSE:240            
##  TRUE :15                             TRUE :6              
##  hospitalizedCumulative  negative       negativeIncrease onVentilatorCumulative
##  Mode :logical          Mode :logical   Mode :logical    Mode :logical         
##  FALSE:246              FALSE:246       FALSE:246        FALSE:225             
##                                                          TRUE :21              
##  onVentilatorCurrently  positive       positiveIncrease recovered      
##  Mode :logical         Mode :logical   Mode :logical    Mode :logical  
##  FALSE:232             FALSE:246       FALSE:246        FALSE:232      
##  TRUE :14                                               TRUE :14       
##    states        totalTestResults totalTestResultsIncrease positivity     
##  Mode :logical   Mode :logical    Mode :logical            Mode :logical  
##  FALSE:246       FALSE:246        FALSE:246                FALSE:246      
## 
xdf %>% dplyr::select(!Gathering) %>% is.na() %>% summary()
##     Date          MaxTempF        AvgTempF        MinTempF      
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246       FALSE:246       FALSE:246      
##  MaxDewPointF    AvgDewPointF    MinDewPointF    MaxHumidityPercent
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical     
##  FALSE:246       FALSE:246       FALSE:246       FALSE:246         
##  AvgHumidityPercent MinHumidityPercent MaxWindSpeedMPH AvgWindSpeedMPH
##  Mode :logical      Mode :logical      Mode :logical   Mode :logical  
##  FALSE:246          FALSE:246          FALSE:246       FALSE:246      
##  MinWindSpeedMPH MaxPressureHg   AvgPressureHg   MinPressureHg  
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246       FALSE:246       FALSE:246      
##  PrecipitationIN  Weekday        GatheringBinary
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:246       FALSE:246       FALSE:246

Comparing Positivity Rate vs New Positive Cases

Positivity rate - According to John Hopkins University, the positivity rate is a new measure developed to identify either trending increases in transmission rates, or overall drops in testing rates, which can indicate a shortage of tests. It measures the ratio of diagnostic tests which are actually positive (i.e. # of positive tests / # of total tests).

The idea behind this test is that a high positivity score should give leaders indication that either infection is increasing at an unusual rate, or not enough tests are available. The World Health Organization has recommended nations adopt a rule of thumb to keep populations under approximately 5% positivity rate before considering re-opening or relaxing social restrictions.

The positivity rate improves upon strictly measuring the total number of cases or positive results by better accounting for new tests. As time goes on, more individuals will be re-tested, including high risk individuals such as health care workers, who may receive tests and positive results more frequently than others. The positivity rate accounts for this by measuring the total number of positive cases against tests, which will show less bias than strictly reporting on positive cases, which can be inflated by duplicate tests and non-unique cases. A better way to control this bias is to measure the people testing positive over the people who have been tested, although this data is much more difficult to collect, due to HIPAA and privacy concerns with linking data to individuals.

Sources https://www.jhsph.edu/covid-19/articles/covid-19-testing-understanding-the-percent-positive.html https://www.mprnews.org/story/2020/08/12/behind-the-numbers-what-does-a-covid19-positivity-rate-really-mean

Below is a plot that shows the differences between new positive cases vs the positivity rate. You can start to see the explosive nature of the positive case metric. I determined the ratio of Louisiana’s population compared to the US population by using census.gov population values from Jul 1, 2019.

#census.gov reported populations on Jul 1, 2019: Louisiana / USA
lapopratio <- 4648794 / 328239523

positivity_df <- data.frame(cbind(la[c('date', 'positivity')],  us['positivity']))
names(positivity_df) = c('date','LA','US')

positive_df <- data.frame(cbind(la[c('date', 'positive')],  us['positiveIncrease']))
names(positive_df) = c('date','LA','US')
positive_df$US = positive_df$US * lapopratio


grid.arrange(
  gather(positivity_df, key = Region, value = 'Positivity Rate', -date) %>% 
    ggplot(aes(date, `Positivity Rate`)) + 
    geom_line(aes(color=Region), size=.75) + 
    labs(title='Louisiana vs US Positivity Rate', 
         x = 'Date', 
         y='Positivity Rate'), 

  gather(positive_df, key = Region, value = 'Positive Tests', -date) %>% 
    ggplot(aes(date, `Positive Tests`)) + 
    geom_line(aes(color=Region), size=.75) + 
    labs(title='Louisiana vs US Daily New Positive Tests', 
         x = 'Date', 
         y='Positive Tests', 
         subtitle = '*US numbers scaled to 1.42% based on LA population density'),
  ncol = 2)

Evaluation Methods

To compare the prediction performance of our models to reality, we want to train a model and hold back a test set of a given period of days that I’ll call a horizon. We will sum the squared difference between each predicted value and actual value, and then take the average to generate the Average Squared Error (ASE). This metric will be the base comparison for all models, and generally the model with the lowest ASE score is the best performing model.

Helper functions

To help score the models, I have written a couple of helper functions. A takes a fitted model, predictions and the desired horizon, and it generates plots of actual, predicted and confidence intervals and the ASE for the given model.

eval_model <- function(response, predictions, pred_ul = NA, pred_ll = NA, model_name, AIC_val = 0, ending_point = length(response)) {
  num_predictions = length(predictions)
  test_stop <- length(response)
  test_start <- test_stop - num_predictions + 1
  compare_stop <- test_start - 1
  compare_start <- compare_stop - num_predictions + 1
  ASE <- mean((predictions - response[test_start:test_stop])^2)

  # Build predictions dataframe
  df <- data.frame('Predicted' = predictions)
  df$Day = row(df)
  df$Actual = response[test_start:test_stop]
  df <- gather(df, key='Type', value, -Day)
  
  #if we have enough data to plot num_predictions * 2, do it, else use num_predictions
  starting_point <- compare_start - num_predictions + 1
  plot_start <- ifelse(starting_point < 0, compare_start, starting_point)
  day_multiplier <- ifelse(starting_point < 0, -1,-2)
  
  # Build predicted vs actual dataframe
  df <- rbind(df, 
    data.frame("Day"=c(((num_predictions-1)*day_multiplier):0), 
               "Type" = 'Actual', 
               value = response[plot_start:compare_stop]))

  # Built UL/LL dataframes
  ul <- data.frame("Day"=c(1:num_predictions), pred_ul)
  ll <- data.frame("Day"=c(1:num_predictions), pred_ll)
  
  # Build Plot
  comparison_plot <- ggplot() + 
    geom_line(data=df, aes(Day + ending_point - num_predictions, value, color=Type)) + 
    geom_point(size=.75) + 
    labs(title=paste(model_name, 'Performance Evaluation'),
         subtitle=paste0(num_predictions,'-Day Forecast'), 
         x='Day', 
         y='Positivity Rate',
         caption=paste0('ASE: ',round(ASE,6),
                       '\nAIC: ',round(AIC_val,6)))
  
  # Add confidence intervals if supplied
  if (length(pred_ul) == length(predictions)){
    comparison_plot = comparison_plot + 
      geom_line(aes(ul$Day + ending_point - num_predictions, ul$pred_ul), 
                color='grey70', linetype = "dashed") 
  }
  if (length(pred_ll) == length(predictions)){
    comparison_plot = comparison_plot + 
      geom_line(aes(ll$Day + ending_point - num_predictions, ll$pred_ll), 
                color='grey70', linetype = "dashed") 
  }
  
  return(comparison_plot)
}
########## example: ARIMA(12,2) ########## 
e <- est.arma.wge(us$positivity, p=12, q=2)
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta=e$theta,  n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, preds$ul, preds$ll, 'ARMA(12,2)', AIC_val = e$aic) 

preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta=e$theta,  n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, preds$ul, preds$ll, 'ARMA(12,2)', AIC_val = e$aic) 

Another function calculates the rolling ASE, which is a measure of the ASE across several windows in time, for a fitted model. This function uses the evaluation function to plot each window and capture score, and then stores them in a list. This list is then averaged to find the average ASE across several time periods. By default I use a month of training data to evaluate a seven day period. The function finally returns a plot of all windows and actual values to visualize performance across the range of data.

rolling_ASE <- function (df, fitted_model, d=0, s=0, horizon=7, training_size=30, model_name, model_type = 'ARUMA', p, df_XDF=NA){
  ASE = list(ASE = c(), plots = c(), multiplot = NA)
  comp_df <- df %>% dplyr::select(date, positivity)
  comp_df$preds = NA
  names(comp_df) = c('date','Actual','Predicted')
  test_stop <- length(df$positivity)
  loop_end <- floor(test_stop/(training_size+horizon))
  
  for (x in 1:loop_end){
    test_start <- test_stop - horizon + 1
    train_start <- test_start - training_size 
    train_stop <- test_start - 1
    print(paste0('test window: ',test_start,':',test_stop,
                 ', train window: ',train_start,':',train_stop))
    data_window <- df$positivity[train_start:train_stop]
    
    if(model_type == 'ARUMA') {
      preds <- fore.aruma.wge(data_window, 
                              phi=fitted_model$phi, 
                              theta=fitted_model$theta, 
                              s=s, 
                              d=d, 
                              n.ahead = horizon, 
                              lastn = F, 
                              limits = F)
      pred_object <- preds$f
    }
    if(model_type == 'SigPlusNoise') {
      preds <- fore.sigplusnoise.wge(data_window, max.p = p, n.ahead = horizon, limits=F)
      pred_object <- preds$f
    }
    
    if(model_type == 'NNFOR') {
      ts_la <- ts(data_window, start = '1')
      mlp_model = mlp(ts_la, lags = horizon, hd.auto.type = 'cv')
      ?mlp
      preds <- predict(mlp_model, horizon)
      preds$ul <- NA
      preds$ll <- NA
      pred_object <- as.numeric(preds$mean)
    }
    
    if(model_type == 'VAR') {
      vfit=VAR(cbind(Positivity = df$positivity, df_XDF)[train_start:train_stop,], p=p, type='both', season = s)
      preds=predict(vfit,n.ahead=7)
      pred_object <- preds$fcst$Positivity[,1]
      preds$ul <- preds$fcst$Positivity[,3]
      preds$ll <- preds$fcst$Positivity[,2]
    }

    a <- mean((pred_object - df$positivity[test_start:test_stop])^2)
    print(paste('Window ASE:', a))
    ASE$ASE[x] <- a
    comp_df$Predicted[test_start:test_stop] = pred_object

    ASE$plots[x] <-
      plot(eval_model(
        data_window,
        pred_object, 
        model_name = model_name, 
        AIC_val = ifelse(model_type == 'ARUMA', fitted_model$aic, 0), 
        pred_ul = preds$ul, 
        pred_ll = preds$ll,
        ending_point = test_stop))
    test_stop = test_stop - training_size
    
  }
  
  ASE$multiplot <- plot(gather(comp_df, key = Type, value = 'Positivity_Rate', -date) %>% 
         ggplot(aes(date, Positivity_Rate, color=Type)) + geom_line() +
         labs(
           title=paste(model_name, 'Performance Evaluation'),
           subtitle=paste0(horizon,'-Day Forecast Rolling Window'), 
           x='Day', 
           y='Positivity Rate',
           caption=paste0('Mean ASE: ',round(mean(ASE$ASE),6))
         )
       ) 
  return(ASE)
}
########## example: ARIMA(12,2) ########## 
e <- est.arma.wge(us$positivity, p=12, q=2)
test <- rolling_ASE(us, e, s=12, d=1, horizon=7, model_name = 'ARMA(12,2)')

Univariate Analysis

The steps for building and evaluating our models will proceed as follows: - Plot and examine the data - Realization - Spectral Density - Autocorrelation Function - Evaluate stationarity assumptions. There are three stationary conditions: - Mean does not depend on time - Variance is finite and does not depend on time - Correlation only depends on how far apart observations are in time, not where they are - Build models - Check that serial correlation has been removed (AIC scores, ACF, PACF plots) - Check residuals for white noise - Evaluate model forecast performance

Louisiana

Plot the Data:

# Plot the data
plotts.sample.wge(la$positivity)

## $autplt
##  [1] 1.0000000 0.9393071 0.8796384 0.8605446 0.8399214 0.8267488 0.8490105
##  [8] 0.8610586 0.8121303 0.7578559 0.7277727 0.6986874 0.6752319 0.6832333
## [15] 0.6821986 0.6292681 0.5668868 0.5319428 0.4966430 0.4664633 0.4593505
## [22] 0.4569352 0.4031459 0.3462679 0.3230427 0.2910256
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  13.9602603  17.5611464  11.9599796   7.6668080   6.4697337   1.6128896
##   [7]  -2.8203865  -5.0586421  -9.6952719  -7.3153505 -23.1165770  -8.5853171
##  [13]  -6.4269213 -10.4512454  -3.9370940  -6.6895400  -5.0361316  -8.8174379
##  [19]  -8.8408129  -7.7524816 -12.7800934 -10.8831615 -12.3443173  -7.4187515
##  [25] -13.2201935 -12.5483706 -12.6072411 -10.4639780 -10.3061541  -5.7896562
##  [31]  -6.9527205  -4.5729994  -2.6565916 -11.1992997   4.9222172  -3.5726203
##  [37] -17.5737307  -6.4938454  -6.9104502 -10.1859592  -7.9690657 -14.6326883
##  [43] -12.8126497 -12.2660625  -9.3597016 -18.1927808 -10.9197193 -15.3687381
##  [49] -12.2455703 -14.3938321 -14.3987676  -8.0215954 -11.2276907 -13.8151299
##  [55] -18.4209672 -14.1462588 -13.9178056 -10.3542921 -18.1090269  -9.4175468
##  [61] -11.6144991 -12.6302015 -16.0793212 -17.8525413 -14.4047311 -13.6480654
##  [67] -19.4031903  -6.2595560  -6.4716425  -0.1702762  -7.4453558 -16.3967146
##  [73] -11.3842620 -12.6276061 -23.1800642 -15.0448944 -12.1965307 -17.0067301
##  [79] -13.2831602 -15.2489091 -11.9935760 -10.7178034 -13.0921498 -25.3789882
##  [85] -14.9455352 -23.1102712 -24.1520012 -20.1648753 -14.8792309 -16.2779671
##  [91] -15.7492788 -11.4820211 -17.6783892 -34.5124466 -17.2522571 -14.1063598
##  [97] -22.9668371 -28.1216282 -22.5894863 -17.0841672 -18.0459368 -11.7319206
## [103] -14.1578096 -16.1153463 -14.0821580 -12.1054220 -21.5900996 -29.1884372
## [109] -31.9989494 -22.2731413 -16.8603667 -42.1727225 -30.0127617 -10.3486417
## [115] -15.6342856 -16.8323448 -23.1519542 -18.4824601 -20.6476432 -22.7906872
## [121] -16.9156395 -31.3160094 -23.1540488
## 
## $dbz
##   [1]  12.566881  12.293334  11.835242  11.189378  10.351358   9.315923
##   [7]   8.077572   6.631985   4.979074   3.129376   1.116753  -0.979425
##  [13]  -3.007511  -4.747472  -6.005907  -6.756174  -7.136643  -7.314196
##  [19]  -7.399551  -7.449649  -7.491625  -7.533833  -7.564111  -7.546353
##  [25]  -7.426789  -7.156016  -6.718379  -6.146017  -5.505068  -4.868567
##  [31]  -4.297108  -3.833194  -3.503840  -3.325433  -3.307926  -3.457655
##  [37]  -3.778856  -4.274127  -4.943857  -5.784457  -6.784834  -7.920367
##  [43]  -9.144086 -10.377572 -11.511317 -12.430793 -13.068592 -13.440575
##  [49] -13.624521 -13.706896 -13.748785 -13.782319 -13.821469 -13.872769
##  [55] -13.939776 -14.020358 -14.098741 -14.136324 -14.068719 -13.820057
##  [61] -13.338769 -12.633433 -11.773842 -10.856577  -9.969285  -9.175272
##  [67]  -8.514349  -8.009317  -7.672104  -7.507927  -7.517585  -7.698307
##  [73]  -8.043468  -8.541337  -9.172968  -9.909643 -10.711031 -11.526432
##  [79] -12.302013 -12.994333 -13.583945 -14.078851 -14.503984 -14.884222
##  [85] -15.231765 -15.543263 -15.806171 -16.010199 -16.157419 -16.264982
##  [91] -16.359111 -16.464519 -16.594649 -16.745743 -16.895822 -17.009646
##  [97] -17.050528 -16.996280 -16.851085 -16.645161 -16.423033 -16.229091
## [103] -16.097778 -16.049668 -16.091173 -16.215626 -16.404893 -16.631990
## [109] -16.865795 -17.078322 -17.253005 -17.390290 -17.507267 -17.631541
## [115] -17.792797 -18.015496 -18.313830 -18.688042 -19.120778 -19.573264
## [121] -19.984024 -20.276769 -20.383632
# Look deeper at spectral density
parzen.wge(la$positivity, trun=100)
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $pzgram
##   [1]  14.8637790  14.2011980  12.6485602  10.0671892   6.8101322   3.6227123
##   [7]   0.7036374  -2.0633900  -4.3842228  -6.5433414  -8.2611328  -8.4389751
##  [13]  -7.6413777  -6.8009486  -6.2911189  -6.1602063  -6.3694419  -6.8987848
##  [19]  -7.7062083  -8.6863538  -9.5077432  -9.6557081  -9.2342817  -8.9830525
##  [25]  -9.2278061  -9.6798562  -9.7913461  -9.1393431  -7.8987016  -6.5685260
##  [31]  -5.5392348  -4.6515485  -2.9494978  -0.8164328   0.2786933  -0.3491239
##  [37]  -2.7974368  -6.3739983  -8.8812547 -10.0296933 -11.2356596 -12.5004769
##  [43] -13.5470125 -13.8931208 -13.8616074 -14.0405632 -14.3322121 -14.6024260
##  [49] -14.6257072 -13.9303719 -12.8233618 -12.2324784 -12.6920818 -14.1670234
##  [55] -15.7658464 -15.9690501 -15.1426719 -14.4520904 -14.0367957 -13.7506921
##  [61] -13.7890897 -14.4329786 -15.5991766 -16.6962898 -16.3502735 -14.2029150
##  [67] -11.2635191  -7.8543781  -5.0149509  -3.6453479  -4.0326645  -6.1283961
##  [73]  -9.2993396 -11.9398547 -13.3944641 -13.8986525 -13.5298106 -13.2516942
##  [79] -13.2951733 -13.1402295 -12.7646338 -12.8686582 -13.8904306 -15.5765112
##  [85] -17.3015466 -18.7823532 -19.4344998 -18.5676263 -16.9718976 -15.3765724
##  [91] -14.3372998 -14.4082388 -15.7083639 -17.3837409 -18.0321588 -18.3635612
##  [97] -19.4736057 -20.5164707 -19.8122370 -18.0227462 -16.4113405 -15.4494854
## [103] -14.9042987 -14.1846154 -13.5910594 -14.0910964 -16.3214414 -20.2716849
## [109] -23.0386491 -22.0845605 -20.6856531 -18.5019235 -16.0286579 -14.7041570
## [115] -15.0159308 -16.8615493 -19.4498201 -21.0946956 -21.2866714 -21.0455908
## [121] -21.2006319 -22.0250180 -22.6403932

There is strong evidence of serial correlation in this data as can be seen in the peaks in the spectral density plots and the slowly damping, sinusoidal behavior in the autocorrelation function plots. There is significant evidence of a weekly trend, as can be seen in the Parzen Window (spectral density plot). The frequency peak at .14 indicates a weekly trend (1/7 = .143), with potentially an echo at .28. Increasing the truncation point of the Parzen function helps highlight some additional points that suggest possibly a monthly factor as well.

Evaluate Stationarity:

There are several indications that suggest the data is NOT stationary, the most obvious being the strong weekly trend with strong peaks each weekend, which effectively increases the mean on weekends. There is a difference in variance which can be seen in the early weeks of the pandemic that have tapered off as time progresses.

Model Construction:

I start by differencing the data, which filters long-term trending behavior and leaves cyclical patterns. Differencing once leaves a good bit of cyclical pattern in the ACF, especially at lag 7.

Identify correlation structures:

Here we are going to look at different filters to try to model out some of the correlation in the data. We’ll also look at a factor table to see if we can identify a weekly or monthly trend.

The overfit table shows strong evidence of 1-B, 1+.445B+B^2, 1-1.247B+B^2 and 1+.802B+B^2 being present in the model. This is more than enough evidence to proceed with the weekly seasonal model.

# Take out 1-B
la.d1 <- artrans.wge(la$positivity,1)

# Taking out another 1-B.. looks like it isn't enough
artrans.wge(la.d1,1)

##   [1] -0.12001810013  0.11197961593  0.01148788128 -0.09768340278  0.08729785421
##   [6] -0.01424222648 -0.00245022335  0.02856748199  0.03663734739 -0.03020620067
##  [11] -0.08490402636  0.07586325598 -0.04411939164  0.01405163464  0.00603443239
##  [16]  0.05326455229 -0.07040279424 -0.01703609861 -0.00895351037 -0.00138728940
##  [21]  0.03352680899 -0.02059646162  0.08058041871 -0.01655174992 -0.12569990013
##  [26]  0.07715699498  0.00970828146 -0.02828671007  0.04922032309 -0.00137427599
##  [31]  0.01757626942 -0.08886926099  0.04368065913 -0.01926462275  0.01546202308
##  [36]  0.01132636604  0.01489236125  0.02424383991 -0.09006170447  0.04445552513
##  [41]  0.01603133393 -0.01939040709  0.00042712593  0.01826882185 -0.00377191437
##  [46] -0.01780837066 -0.00002107674  0.02992415119 -0.03931100761  0.02461497660
##  [51]  0.01422652210 -0.01024811730 -0.03048764593  0.01974597812  0.00900167078
##  [56]  0.00819385113 -0.03064462453  0.00950030564  0.04636873235 -0.05359735437
##  [61]  0.00719330921  0.00019993853  0.02416632150 -0.01474365246 -0.01262301520
##  [66]  0.04615993257 -0.05738147393  0.00585647802  0.02531654811 -0.00093831949
##  [71] -0.01817074820  0.04329415542 -0.04213960697  0.00626173920 -0.00566945709
##  [76]  0.00745556817  0.00195300117  0.00524211825  0.00262408214  0.00460223860
##  [81] -0.02713178495  0.00522217386  0.01741746191 -0.01115368628  0.00752134400
##  [86]  0.00774049302  0.01478958929 -0.06086352002  0.03782370394 -0.00571849966
##  [91]  0.00319051199  0.00398135599  0.04963294085 -0.05888886023 -0.03562008158
##  [96]  0.04411647176  0.00224009636  0.00153450585 -0.00789657627  0.03008822968
## [101] -0.01229604121 -0.05886386305  0.03469231074  0.00710503644  0.00257321289
## [106]  0.00022733426  0.04154065525 -0.04786083785 -0.01709010640  0.01517072373
## [111]  0.00760007924  0.00649024584  0.02373510913 -0.00284543631 -0.01272601102
## [116] -0.08239005417  0.04713628416  0.02718539078 -0.00738986440 -0.00059194356
## [121]  0.03070118614 -0.02070185623 -0.05368870817  0.02827825825  0.01409853916
## [126]  0.00576960680 -0.00084032111  0.02024613738 -0.01441860235 -0.03524539380
## [131]  0.01209539642  0.01070138410  0.00084646748  0.00828549360  0.01419305283
## [136] -0.00792484179 -0.04416309640  0.00375735815  0.03480591454 -0.00913183129
## [141]  0.00713310866  0.01075547427 -0.01421246068 -0.03455009149  0.02639598188
## [146]  0.00021892746  0.01219897879 -0.01321799746  0.03121307054 -0.01499341370
## [151] -0.04776882472  0.02695663191  0.00674473559  0.01212144857 -0.00504886828
## [156]  0.03724066022 -0.06437056017 -0.00005126073  0.02327608778  0.00336498625
## [161] -0.00103675954  0.00310067789  0.01971938490 -0.02913216776 -0.01216822883
## [166]  0.01335914467 -0.00013283852  0.00748879028  0.01075283531 -0.01444866773
## [171]  0.01637641523 -0.04955039165  0.02621505529  0.01434378900 -0.01189553928
## [176]  0.00512143710  0.01253283708 -0.00618458919 -0.03496356980  0.03070706559
## [181] -0.00885731849  0.00400893776  0.00250957556  0.03818324918 -0.03313117378
## [186] -0.04599125522  0.04620206871 -0.00982100992  0.00260166582  0.00312275101
## [191]  0.02014836905 -0.02033125343 -0.02344872791  0.01835402238  0.01194806885
## [196] -0.00898995412  0.00283657842  0.02190296124 -0.02927985203 -0.01345968115
## [201]  0.01367358904  0.00621727754  0.00114851488  0.00516535494  0.01112481513
## [206] -0.01896296615 -0.02146644421  0.01711274737  0.01075193799 -0.00423858946
## [211] -0.00013615570  0.04181976478 -0.04444093558 -0.03623206054  0.03222659037
## [216]  0.01219654120 -0.00589440433  0.00152342606  0.02251320599 -0.02345323556
## [221] -0.03253380372  0.02430159488  0.01522738477 -0.00666540167  0.00040987785
## [226]  0.01831575860 -0.01954011115 -0.02892179674  0.03014970836  0.00318487094
## [231]  0.00191084383 -0.00973136741  0.03725493099 -0.03444842445 -0.04309372782
## [236]  0.04491734584  0.01680699336 -0.01430450447 -0.00035522384  0.02947585241
## [241] -0.00995478282 -0.07682562404  0.04821360473  0.01244096381
# strong evidence of a 7 day cycle + possibly a monthly cycle
la.s7 <- artrans.wge(la$positivity, c(rep(0,6),1))

# Add (1-B) + s7
la.d1.s7 <- artrans.wge(la.s7,1)

# Try to model weekly + monthly data
la.s7.s12 <- artrans.wge(la$positivity, c(rep(0,6),1,0,0,0,0,1))

# Overfit table
est.ar.wge(la$positivity, p=15)
## 
## Coefficients of Original polynomial:  
## 0.7229 0.1790 0.2398 -0.1301 -0.0733 0.0620 0.3964 -0.3014 -0.2877 -0.1168 0.0513 0.0441 0.0674 0.3147 -0.1753 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.9959B              1.0041               0.9959       0.0000
## 1-1.2348B+0.9762B^2    0.6324+-0.7902i      0.9881       0.1426
## 1+0.4517B+0.9649B^2   -0.2340+-0.9908i      0.9823       0.2869
## 1-1.8259B+0.8957B^2    1.0193+-0.2785i      0.9464       0.0424
## 1+1.6899B+0.8887B^2   -0.9508+-0.4704i      0.9427       0.4269
## 1+1.1264B+0.7953B^2   -0.7082+-0.8694i      0.8918       0.3588
## 1+0.8673B             -1.1530               0.8673       0.5000
## 1-0.2981B+0.6757B^2    0.2206+-1.1964i      0.8220       0.2210
## 1-0.5036B              1.9857               0.5036       0.0000
##   
##   
## 
## Coefficients of Original polynomial:  
## 0.7229 0.1790 0.2398 -0.1301 -0.0733 0.0620 0.3964 -0.3014 -0.2877 -0.1168 0.0513 0.0441 0.0674 0.3147 -0.1753 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.9959B              1.0041               0.9959       0.0000
## 1-1.2348B+0.9762B^2    0.6324+-0.7902i      0.9881       0.1426
## 1+0.4517B+0.9649B^2   -0.2340+-0.9908i      0.9823       0.2869
## 1-1.8259B+0.8957B^2    1.0193+-0.2785i      0.9464       0.0424
## 1+1.6899B+0.8887B^2   -0.9508+-0.4704i      0.9427       0.4269
## 1+1.1264B+0.7953B^2   -0.7082+-0.8694i      0.8918       0.3588
## 1+0.8673B             -1.1530               0.8673       0.5000
## 1-0.2981B+0.6757B^2    0.2206+-1.1964i      0.8220       0.2210
## 1-0.5036B              1.9857               0.5036       0.0000
##   
## 
## $phi
##  [1]  0.72293214  0.17904609  0.23984475 -0.13008877 -0.07327756  0.06198856
##  [7]  0.39641995 -0.30136793 -0.28766989 -0.11683003  0.05126236  0.04414203
## [13]  0.06741100  0.31465706 -0.17526686
## 
## $res
##   [1] -0.03093945037  0.02316236956 -0.05435390974 -0.00502183438  0.00893904094
##   [6] -0.01985489760  0.02518009493 -0.00323618159 -0.01080959161  0.05450484314
##  [11]  0.05434986643  0.01894490423 -0.02611320488  0.00463390908 -0.02988421379
##  [16] -0.01595229496  0.00519612379  0.02334120053 -0.01637060949  0.02673726934
##  [21] -0.02746859385 -0.02094560817  0.00948456862 -0.02138901615  0.01445368950
##  [26]  0.03925759433 -0.03524684298 -0.01054355581  0.00942017257 -0.01413700563
##  [31]  0.00694955701 -0.04274641357  0.00074863942 -0.01905898380  0.00565668039
##  [36] -0.02889700609 -0.01189626903 -0.01409257323 -0.02753259822  0.00034614558
##  [41] -0.00850118185  0.00066293073  0.01340145691  0.00903591305 -0.02212165781
##  [46] -0.02874659746 -0.03762965113  0.00770188222 -0.00551425507  0.01429056095
##  [51] -0.01787291268 -0.00059923134 -0.00156357314 -0.00991733273 -0.01161547248
##  [56] -0.00896680848 -0.01191706546  0.01956372890 -0.01415446599 -0.02954068870
##  [61]  0.01749759469  0.00517144936  0.00130098304 -0.02678301296  0.00678500035
##  [66]  0.00814910002 -0.01142255741  0.00807402729 -0.01465019447 -0.01094130636
##  [71]  0.00186169458 -0.00131120746 -0.00604921370  0.03148190920 -0.02682804486
##  [76]  0.00412801364 -0.00174523049 -0.00123813447 -0.01123950065  0.00317430085
##  [81] -0.00282394477  0.00263284700 -0.00311323827 -0.00067644099  0.00368140353
##  [86] -0.00288795943  0.00256798041 -0.00542778472  0.02242509916 -0.02085880014
##  [91]  0.00566567188 -0.00770799376  0.00408359714  0.00116957504  0.04409157264
##  [96] -0.00638623605 -0.02161682831 -0.00287760560  0.00660384509  0.01811151176
## [101] -0.00261361664 -0.00116705444  0.01225968475  0.00024633298 -0.00248237358
## [106] -0.00568342337  0.00470324489  0.00321782831  0.01666260609 -0.00190960440
## [111]  0.01333832066  0.00304272435  0.00417870431  0.00846533013  0.02955830254
## [116]  0.00731876833  0.01483670557 -0.03764748296 -0.01929639747  0.00351744585
## [121]  0.01039928781 -0.00854192500  0.00054005230  0.01490819192  0.00938605645
## [126]  0.00770159066 -0.00417226682  0.00046251487 -0.00964401877 -0.00644861392
## [131] -0.00429676443  0.01314316916  0.00481705246 -0.00584066542 -0.00678062295
## [136] -0.00184183800 -0.00395409583  0.00431774407 -0.00443480516 -0.01568071375
## [141]  0.00408165357  0.00177380610  0.00267825774 -0.00603350603 -0.00899924471
## [146] -0.01134564320  0.00916042653 -0.00580768182  0.00418374479 -0.01636468317
## [151]  0.00200478699  0.00222983810 -0.01240003293 -0.00612019352 -0.01237420014
## [156]  0.00393755026 -0.00059175998  0.01909125446 -0.02822531155 -0.00629549502
## [161]  0.00136862439  0.00741735403 -0.00360759578 -0.00595179654 -0.00683535968
## [166] -0.00068930161  0.00915453494  0.00076502854 -0.00714786447 -0.00200711739
## [171]  0.00654471037 -0.02202500232  0.01806906077 -0.01393317717 -0.00412223383
## [176]  0.00358651115 -0.00263152438 -0.00939463150 -0.00261176383  0.00370279100
## [181] -0.00694669263  0.00588280941 -0.00983992258 -0.00498181553 -0.00977576459
## [186]  0.02602126906  0.00301226969 -0.01863472947 -0.00355022985 -0.00860476524
## [191]  0.00329891257 -0.00487912168 -0.00550424536 -0.00582125658  0.00506704835
## [196] -0.00146964615  0.00898300503  0.00374046125 -0.00121913464 -0.00277754085
## [201] -0.00943344986  0.00076418798 -0.00672419290 -0.00331745977  0.00225067615
## [206]  0.00446921729  0.00164313697  0.00131682864 -0.00458535259 -0.00253809423
## [211]  0.00121465831  0.00211380925 -0.00523795138  0.02301279532  0.00714116899
## [216] -0.01666709162 -0.00992813954  0.00061765639  0.00485997939 -0.00195597126
## [221] -0.00276673410  0.00177202064 -0.00053047284  0.00152742254  0.00312650034
## [226]  0.00371745232 -0.00117132795 -0.00724414542 -0.00125482227 -0.00101012899
## [231]  0.00733384230 -0.00006790254  0.00697149684 -0.00651165092  0.01388212638
## [236]  0.00286561357 -0.01701313400 -0.00316294895  0.01119596746  0.01067369287
## [241]  0.00460069718  0.01018875659  0.02827800653 -0.01218119477 -0.00701711434
## [246] -0.01030342246
## 
## $avar
## [1] 0.0002103418
## 
## $aic
## [1] -8.336695
## 
## $aicc
## [1] -7.317654
## 
## $bic
## [1] -8.108706
# Generate a factor table for a (1-B^7)
tswge::factor.wge(c(rep(0,6),1))
## 
## Coefficients of Original polynomial:  
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.0000B              1.0000               1.0000       0.0000
## 1+0.4450B+1.0000B^2   -0.2225+-0.9749i      1.0000       0.2857
## 1-1.2470B+1.0000B^2    0.6235+-0.7818i      1.0000       0.1429
## 1+1.8019B+1.0000B^2   -0.9010+-0.4339i      1.0000       0.4286
##   
## 

Estimate parameters and check residuals:

TSWGE offers an aic5 function that iterates through values of P (phi) and Q (theta), fits a model, and returns the top 5 performing models by either AIC or BIC value. The AIC and BIC penalize models for the number of parameters used, with BIC favoring simpler models even more so than AIC. This gives us a pretty good idea of the best fitting model as a starting point.

The ACF and PACF plots look like white noise, as expected. The Ljung-Box test fails to reject the null hyphothesis that the residuals are white noise. The generated spectral densities look like they have a fairly good fit, the final peak is overpronounced in our generations compared to the actual data. The ACF comparisons are all over the board, but with peaks in the same general spots (lags at 7, 14, 21). The generated realizations look quite a bit different, but the spectral densities and ACF plots look fairly similar, indicating this could be a useful model.

# Estimate the model params
aic5.wge(la.s7, type = 'bic', p=0:20,q=0:2) #15
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Error in aic calculation at 18 1 
## Five Smallest Values of  bic
##       p    q        bic
## 46   15    0  -8.256251
## 49   16    0  -8.249310
## 48   15    2  -8.247515
## 42   13    2  -8.246550
## 52   17    0  -8.242890
e <- est.arma.wge(la.s7, p=15, q=0)
## 
## Coefficients of Original polynomial:  
## 0.6819 0.2329 0.1795 -0.1438 0.0353 0.0725 -0.6438 0.3639 0.0825 -0.0664 -0.1285 0.0148 0.2163 -0.3415 0.3063 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.8995B+0.9868B^2    0.9625+-0.2950i      0.9934       0.0473
## 1+0.9661B+0.9242B^2   -0.5226+-0.8994i      0.9614       0.3338
## 1+1.3916B+0.9115B^2   -0.7634+-0.7172i      0.9547       0.3800
## 1+1.8666B+0.8970B^2   -1.0405+-0.1795i      0.9471       0.4728
## 1-0.9455B              1.0576               0.9455       0.0000
## 1-0.1935B+0.8186B^2    0.1182+-1.0989i      0.9048       0.2329
## 1-1.3104B+0.7391B^2    0.8865+-0.7531i      0.8597       0.1121
## 1-0.5571B+0.7181B^2    0.3879+-1.1145i      0.8474       0.1967
##   
## 
# Ljung Box Test shows white noise residuals
ljung.wge(artrans.wge(la.s7, phi.tr = e$phi))
## Obs 0.03251153 0.005953455 -0.04811914 0.1123665 0.05091219 -0.006740414 0.02838045 0.05305956 -0.01653444 0.03643441 0.01285373 0.02944821 0.008132136 -0.003795011 0.08950678 -0.02009765 0.05161991 0.02098249 0.0549419 0.02685266 -0.081336 0.07919188 -0.08543065 0.0593491
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 24
## 
## $chi.square
## [1] 15.45967
## 
## $df
## [1] 24
## 
## $pval
## [1] 0.9065855
ljung.wge(artrans.wge(la.s7, phi.tr = e$phi), K = 48)

## Obs 0.03251153 0.005953455 -0.04811914 0.1123665 0.05091219 -0.006740414 0.02838045 0.05305956 -0.01653444 0.03643441 0.01285373 0.02944821 0.008132136 -0.003795011 0.08950678 -0.02009765 0.05161991 0.02098249 0.0549419 0.02685266 -0.081336 0.07919188 -0.08543065 0.0593491 -0.03150929 0.05970946 -0.09874056 0.04136814 0.06256806 -0.0520347 -0.04590716 0.06908596 0.04504184 -0.05018569 -0.09408055 -0.05452445 0.1304085 -0.04294655 -0.03968746 -0.02780619 0.0549043 -0.141726 0.1096732 -0.03633432 -0.02496781 -0.03631315 0.05437685 -0.02354808
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 48
## 
## $chi.square
## [1] 44.94805
## 
## $df
## [1] 48
## 
## $pval
## [1] 0.5986569
acf(e$res)
pacf(e$res)
dev.off()
## null device 
##           1
#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(la$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)

for( i in 1: sims)
{
   SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 7, phi = e$phi, plot ="FALSE"), plot = "FALSE")
   lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}


#Compare ACFs
sims = 5
ACF = acf(la$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)

for( i in 1: sims)
{
   ACF2 = acf(gen.aruma.wge(246, s = 7, phi = e$phi, plot = "FALSE"), plot = "FALSE")
   lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}

#Compare Generated Realizations 
eGen = gen.aruma.wge(246, s = 7, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)
## $autplt
##  [1] 1.0000000 0.9736991 0.9626871 0.9439085 0.9248900 0.9046038 0.8813618
##  [8] 0.8742908 0.8401772 0.8238030 0.8047006 0.7929718 0.7858042 0.7776364
## [15] 0.7873074 0.7724358 0.7735335 0.7716518 0.7714489 0.7694762 0.7618503
## [22] 0.7669221 0.7428109 0.7309969 0.7125539 0.6936437
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  18.959905  11.822348   8.431301   5.279660   2.998245   1.920224
##   [7]   2.019772   0.889596  -4.071676   2.360229   3.035087   5.596055
##  [13] -14.718606  -5.706822  -5.895885  -9.281965 -10.965889  -8.360391
##  [19] -11.536625  -8.941074 -11.201632  -9.397568 -10.128222 -11.581647
##  [25] -10.530920 -10.415388 -12.223071 -13.714795 -15.929249 -11.416607
##  [31] -14.650649 -13.000973  -8.652443  -6.858406  -9.184427 -13.538787
##  [37] -13.512267 -13.838175 -11.400413 -13.597119 -16.074031 -17.743589
##  [43] -19.187122 -15.621058 -16.246249 -19.156378 -20.689526 -21.787558
##  [49] -21.207093 -23.609956 -17.776729 -19.532000 -15.296678 -16.045822
##  [55] -15.673240 -20.973239 -18.562394 -17.623196 -16.395773 -16.102934
##  [61] -13.468320 -18.070707 -17.961197 -22.573351 -15.967344 -17.201750
##  [67] -16.657263 -11.192350  -9.885698 -14.675763  -8.173531 -11.008147
##  [73] -18.271583 -31.978990 -23.808425 -24.623456 -29.502731 -24.997538
##  [79] -17.720580 -23.963103 -14.550412 -15.840960 -17.610135 -26.105177
##  [85] -24.648677 -19.517523 -17.956652 -21.135251 -25.198270 -20.907685
##  [91] -23.646234 -19.497703 -25.984641 -20.813359 -28.070662 -21.818354
##  [97] -29.093402 -18.943961 -33.424419 -24.213975 -21.867636 -25.267836
## [103] -20.129915 -18.508756  -3.395199 -15.594638  -9.780505 -12.382323
## [109] -23.644788 -23.055684 -20.305960 -17.995426 -22.721691 -22.945210
## [115] -25.848928 -15.746649 -15.911528 -19.473489 -26.933756 -22.765069
## [121] -17.658643 -23.506246 -22.314047
## 
## $dbz
##   [1]  12.98338914  12.67366873  12.15552791  11.42660225  10.48495692
##   [6]   9.33124119   7.97286582   6.43161294   4.75581155   3.03426865
##  [11]   1.39679121  -0.02681045  -1.18680272  -2.15349125  -3.06123556
##  [16]  -4.02262187  -5.08250916  -6.20823593  -7.30055818  -8.23754381
##  [21]  -8.95433325  -9.49493285  -9.97289978 -10.48822020 -11.07456247
##  [26] -11.68590510 -12.21180282 -12.53019870 -12.58691479 -12.43198324
##  [31] -12.17306092 -11.90820127 -11.69804179 -11.57226267 -11.54552778
##  [36] -11.62965214 -11.83896601 -12.19009842 -12.69849368 -13.37343017
##  [41] -14.21210568 -15.19254148 -16.26576559 -17.35119214 -18.34460413
##  [46] -19.14695255 -19.70176485 -20.00808893 -20.09903760 -20.01672232
##  [51] -19.80703194 -19.52417350 -19.22771230 -18.97079596 -18.78820767
##  [56] -18.68937058 -18.65576573 -18.64168561 -18.58055500 -18.40142305
##  [61] -18.05518080 -17.53712130 -16.88829227 -16.17555161 -15.46788247
##  [66] -14.82260989 -14.28210319 -13.87578779 -13.62336632 -13.53747101
##  [71] -13.62529072 -13.88913125 -14.32585590 -14.92501128 -15.66537773
##  [76] -16.51005025 -17.40168428 -18.26289651 -19.00989574 -19.58185503
##  [81] -19.96951744 -20.21667709 -20.39246380 -20.55966131 -20.75796337
##  [86] -21.00171502 -21.28492705 -21.58899094 -21.89083122 -22.16819344
##  [91] -22.39672151 -22.53526310 -22.50600964 -22.19544773 -21.50929252
##  [96] -20.45624686 -19.16111365 -17.78777125 -16.46736359 -15.28023881
## [101] -14.26794488 -13.44900781 -12.82996291 -12.41168250 -12.19267076
## [106] -12.17060484 -12.34282141 -12.70602136 -13.25515670 -13.98119877
## [111] -14.86726515 -15.88261404 -16.97503703 -18.06562416 -19.05583463
## [116] -19.85709990 -20.43230235 -20.81111603 -21.06017500 -21.23853777
## [121] -21.37384042 -21.46458579 -21.49722005
plotts.sample.wge(la$positivity)
## $autplt
##  [1] 1.0000000 0.9393071 0.8796384 0.8605446 0.8399214 0.8267488 0.8490105
##  [8] 0.8610586 0.8121303 0.7578559 0.7277727 0.6986874 0.6752319 0.6832333
## [15] 0.6821986 0.6292681 0.5668868 0.5319428 0.4966430 0.4664633 0.4593505
## [22] 0.4569352 0.4031459 0.3462679 0.3230427 0.2910256
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  13.9602603  17.5611464  11.9599796   7.6668080   6.4697337   1.6128896
##   [7]  -2.8203865  -5.0586421  -9.6952719  -7.3153505 -23.1165770  -8.5853171
##  [13]  -6.4269213 -10.4512454  -3.9370940  -6.6895400  -5.0361316  -8.8174379
##  [19]  -8.8408129  -7.7524816 -12.7800934 -10.8831615 -12.3443173  -7.4187515
##  [25] -13.2201935 -12.5483706 -12.6072411 -10.4639780 -10.3061541  -5.7896562
##  [31]  -6.9527205  -4.5729994  -2.6565916 -11.1992997   4.9222172  -3.5726203
##  [37] -17.5737307  -6.4938454  -6.9104502 -10.1859592  -7.9690657 -14.6326883
##  [43] -12.8126497 -12.2660625  -9.3597016 -18.1927808 -10.9197193 -15.3687381
##  [49] -12.2455703 -14.3938321 -14.3987676  -8.0215954 -11.2276907 -13.8151299
##  [55] -18.4209672 -14.1462588 -13.9178056 -10.3542921 -18.1090269  -9.4175468
##  [61] -11.6144991 -12.6302015 -16.0793212 -17.8525413 -14.4047311 -13.6480654
##  [67] -19.4031903  -6.2595560  -6.4716425  -0.1702762  -7.4453558 -16.3967146
##  [73] -11.3842620 -12.6276061 -23.1800642 -15.0448944 -12.1965307 -17.0067301
##  [79] -13.2831602 -15.2489091 -11.9935760 -10.7178034 -13.0921498 -25.3789882
##  [85] -14.9455352 -23.1102712 -24.1520012 -20.1648753 -14.8792309 -16.2779671
##  [91] -15.7492788 -11.4820211 -17.6783892 -34.5124466 -17.2522571 -14.1063598
##  [97] -22.9668371 -28.1216282 -22.5894863 -17.0841672 -18.0459368 -11.7319206
## [103] -14.1578096 -16.1153463 -14.0821580 -12.1054220 -21.5900996 -29.1884372
## [109] -31.9989494 -22.2731413 -16.8603667 -42.1727225 -30.0127617 -10.3486417
## [115] -15.6342856 -16.8323448 -23.1519542 -18.4824601 -20.6476432 -22.7906872
## [121] -16.9156395 -31.3160094 -23.1540488
## 
## $dbz
##   [1]  12.566881  12.293334  11.835242  11.189378  10.351358   9.315923
##   [7]   8.077572   6.631985   4.979074   3.129376   1.116753  -0.979425
##  [13]  -3.007511  -4.747472  -6.005907  -6.756174  -7.136643  -7.314196
##  [19]  -7.399551  -7.449649  -7.491625  -7.533833  -7.564111  -7.546353
##  [25]  -7.426789  -7.156016  -6.718379  -6.146017  -5.505068  -4.868567
##  [31]  -4.297108  -3.833194  -3.503840  -3.325433  -3.307926  -3.457655
##  [37]  -3.778856  -4.274127  -4.943857  -5.784457  -6.784834  -7.920367
##  [43]  -9.144086 -10.377572 -11.511317 -12.430793 -13.068592 -13.440575
##  [49] -13.624521 -13.706896 -13.748785 -13.782319 -13.821469 -13.872769
##  [55] -13.939776 -14.020358 -14.098741 -14.136324 -14.068719 -13.820057
##  [61] -13.338769 -12.633433 -11.773842 -10.856577  -9.969285  -9.175272
##  [67]  -8.514349  -8.009317  -7.672104  -7.507927  -7.517585  -7.698307
##  [73]  -8.043468  -8.541337  -9.172968  -9.909643 -10.711031 -11.526432
##  [79] -12.302013 -12.994333 -13.583945 -14.078851 -14.503984 -14.884222
##  [85] -15.231765 -15.543263 -15.806171 -16.010199 -16.157419 -16.264982
##  [91] -16.359111 -16.464519 -16.594649 -16.745743 -16.895822 -17.009646
##  [97] -17.050528 -16.996280 -16.851085 -16.645161 -16.423033 -16.229091
## [103] -16.097778 -16.049668 -16.091173 -16.215626 -16.404893 -16.631990
## [109] -16.865795 -17.078322 -17.253005 -17.390290 -17.507267 -17.631541
## [115] -17.792797 -18.015496 -18.313830 -18.688042 -19.120778 -19.573264
## [121] -19.984024 -20.276769 -20.383632
# Check performance
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 7, lastn = T, limits = F)
eval_model(la$positivity,preds$f, model_name = 'AR(15) With Weekly Trend', AIC_val =  e$aic) #ASE .000647

preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 90, lastn = T, limits = F)
eval_model(la$positivity,preds$f, model_name = 'AR(15) With Weekly Trend', AIC_val =  e$aic) #ASE .000202

rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'AR(15) with Weekly Trend') #ASE .000277
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000646772835406658"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.00011534695778759"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000625442246018444"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000474435722736831"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000099788150409411"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000265503096318144"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).

Fitting additional models:

Weekly ARUMA(14,1,0)

ARUMA(14,1,0) with a weekly trend shows a slight improvement in the 7 day prediction, but a very poor performance in the 90 day trend. This is due to the 1-B term induced by the difference, that causes the trend to continue in a downward trend seen in the training data.

aic5.wge(la.d1.s7, p=0:15, q=0:2, type='bic') #BIC recommends 14,0
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 43   14    0  -8.248609
## 46   15    0  -8.234827
## 44   14    1  -8.230948
## 45   14    2  -8.228824
## 47   15    1  -8.214616
e <- est.arma.wge(la.d1.s7,p = 14, q=0)
## 
## Coefficients of Original polynomial:  
## -0.2868 -0.0488 0.1346 -0.0076 0.0328 0.1071 -0.5340 -0.1427 -0.0600 -0.1284 -0.2537 -0.2335 -0.0086 -0.3409 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.8991B+0.9887B^2    0.9604+-0.2985i      0.9944       0.0480
## 1+0.9691B+0.9260B^2   -0.5233+-0.8978i      0.9623       0.3340
## 1+1.3945B+0.9133B^2   -0.7634+-0.7156i      0.9557       0.3801
## 1+1.8695B+0.8998B^2   -1.0388+-0.1793i      0.9486       0.4728
## 1-0.1828B+0.8139B^2    0.1123+-1.1027i      0.9022       0.2339
## 1-1.3169B+0.7598B^2    0.8665+-0.7518i      0.8717       0.1137
## 1-0.5477B+0.7325B^2    0.3738+-1.1070i      0.8558       0.1982
##   
## 
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta=e$theta, d=1, s=7, n.ahead = 7, lastn = T, limits = F)

eval_model(la$positivity,preds$f, preds$ul, preds$ll,
           'ARIMA(14,1,0) With Weekly Trend', AIC_val = e$aic) #ASE = .000541

preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta=e$theta, d=1, s=7, n.ahead = 90, lastn = T, limits = F)

eval_model(la$positivity,preds$f, preds$ul, preds$ll,
           'ARIMA(14,1,0) With Weekly Trend', AIC_val = e$aic) #ASE = .014017  (BAD)

rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'ARIMA(14,1,0) with Weekly Trend') #ASE .000679
## [1] "test window: 240:246, train window: 210:239"

## [1] "Window ASE: 0.00127641370441219"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"

## [1] "Window ASE: 0.000113147657616594"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"

## [1] "Window ASE: 0.000384666162028781"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"

## [1] "Window ASE: 0.00145540483664713"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"

## [1] "Window ASE: 0.000502713903767943"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"

## [1] "Window ASE: 0.000340185810696027"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

Signal Plus Noise

The signal plus noise model performs very well in the 7 day window, but very poorly over the 90 day window. This is a deterministic signal with a stationary mean model, and our results aren’t surprising since it doesn’t account for the weekly trend well.

########## Sig Plus Noise ########## 
preds <- fore.sigplusnoise.wge(la$positivity, max.p = 15, n.ahead = 7, limits=F)
## 
## Coefficients of Original polynomial:  
## 0.7543 0.0962 0.2365 -0.2347 -0.0441 0.1686 0.3945 -0.3197 -0.2312 -0.0513 0.0608 0.0299 0.0172 0.2035 -0.1219 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.2292B+0.9629B^2    0.6383+-0.7944i      0.9813       0.1423
## 1-0.9726B              1.0282               0.9726       0.0000
## 1+0.4746B+0.9414B^2   -0.2521+-0.9993i      0.9703       0.2893
## 1+1.6350B+0.8406B^2   -0.9725+-0.4938i      0.9168       0.4252
## 1-1.6901B+0.7766B^2    1.0882+-0.3218i      0.8812       0.0458
## 1+1.0723B+0.7471B^2   -0.7177+-0.9075i      0.8643       0.3565
## 1+0.8363B             -1.1958               0.8363       0.5000
## 1-0.3085B+0.5926B^2    0.2603+-1.2727i      0.7698       0.2179
## 1-0.5721B              1.7479               0.5721       0.0000
##   
## 

eval_model(la$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise', 0) #ASE = .000385

preds <- fore.sigplusnoise.wge(la$positivity, max.p = 15, n.ahead = 90, limits=F)
## 
## Coefficients of Original polynomial:  
## 0.7576 0.0956 0.2440 -0.2442 -0.0649 0.1687 0.4054 -0.3083 -0.2337 -0.0644 0.0805 0.0514 0.0181 0.1827 -0.1313 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.2254B+0.9632B^2    0.6361+-0.7960i      0.9814       0.1427
## 1-0.9704B              1.0305               0.9704       0.0000
## 1+0.4720B+0.9372B^2   -0.2518+-1.0018i      0.9681       0.2892
## 1+1.6345B+0.8382B^2   -0.9750+-0.4924i      0.9155       0.4256
## 1-1.6494B+0.7501B^2    1.0995+-0.3527i      0.8661       0.0494
## 1+1.0698B+0.7462B^2   -0.7168+-0.9090i      0.8639       0.3563
## 1+0.8339B             -1.1992               0.8339       0.5000
## 1-0.2651B+0.5828B^2    0.2275+-1.2900i      0.7634       0.2222
## 1-0.6575B              1.5209               0.6575       0.0000
##   
## 

eval_model(la$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise', 0) #ASE = .003437 (BAD)

rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'SigPlusNoise', model_type = 'SigPlusNoise', p = 15) #ASE .000418
## [1] "test window: 240:246, train window: 210:239"
## 
## Coefficients of Original polynomial:  
## 0.4613 -0.6072 0.2538 -0.2611 -0.0181 -0.1028 0.3613 -0.1860 -0.1092 -0.0082 -0.2988 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.2186B+0.9706B^2    0.6278+-0.7976i      0.9852       0.1439
## 1+0.4791B+0.9365B^2   -0.2558+-1.0012i      0.9677       0.2898
## 1-1.7268B+0.8266B^2    1.0445+-0.3446i      0.9092       0.0507
## 1+1.4620B+0.8123B^2   -0.8999+-0.6490i      0.9013       0.4006
## 1-0.1951B+0.6631B^2    0.1471+-1.2191i      0.8143       0.2309
## 1+0.7381B             -1.3549               0.7381       0.5000
##   
## 

## [1] "Window ASE: 0.000799097250179986"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"
## 
## Coefficients of Original polynomial:  
## 0.2158 -0.3047 -0.1716 -0.2450 -0.2585 -0.0610 0.5038 -0.3538 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+0.4915B+0.9561B^2   -0.2570+-0.9899i      0.9778       0.2904
## 1-1.1763B+0.9488B^2    0.6199+-0.8183i      0.9741       0.1468
## 1+1.6727B+0.8580B^2   -0.9748+-0.4640i      0.9263       0.4293
## 1-1.2037B+0.4546B^2    1.3239+-0.6686i      0.6742       0.0744
##   
## 

## [1] "Window ASE: 0.000222706814850399"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## 
## Coefficients of Original polynomial:  
## 0.1201 -0.2251 -0.1223 -0.1370 -0.2003 0.2381 0.3289 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.1803B+0.9519B^2    0.6199+-0.8162i      0.9757       0.1466
## 1+0.5098B+0.8123B^2   -0.3138+-1.0642i      0.9013       0.2956
## 1-0.8070B              1.2392               0.8070       0.0000
## 1+1.3574B+0.5270B^2   -1.2878+-0.4888i      0.7260       0.4423
##   
## 

## [1] "Window ASE: 0.000112207413882515"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"
## 
## Coefficients of Original polynomial:  
## 0.6905 -0.4220 0.1930 -0.2575 -0.0349 0.1692 0.3951 -0.3312 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.1862B+0.9627B^2    0.6161+-0.8119i      0.9812       0.1467
## 1+0.4931B+0.9174B^2   -0.2687+-1.0089i      0.9578       0.2914
## 1+1.4515B+0.6808B^2   -1.0660+-0.5766i      0.8251       0.4211
## 1-1.4488B+0.5508B^2    1.3152+-0.2928i      0.7422       0.0349
##   
## 

## [1] "Window ASE: 0.000207049686158031"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"
## 
## Coefficients of Original polynomial:  
## 0.5644 -0.5562 0.0205 -0.1355 -0.0901 -0.0984 0.1800 -0.0990 -0.3832 0.2435 -0.3451 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.1705B+0.9263B^2    0.6318+-0.8249i      0.9624       0.1460
## 1+0.3801B+0.9120B^2   -0.2084+-1.0262i      0.9550       0.2819
## 1+1.4687B+0.8672B^2   -0.8468+-0.6604i      0.9312       0.3946
## 1-1.7165B+0.8443B^2    1.0165+-0.3888i      0.9189       0.0581
## 1+0.8842B             -1.1310               0.8842       0.5000
## 1-0.4104B+0.6311B^2    0.3251+-1.2161i      0.7944       0.2084
##   
## 

## [1] "Window ASE: 0.000603111475107872"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## 
## Coefficients of Original polynomial:  
## -0.1500 -0.0283 0.0601 -0.1501 -0.3244 0.0988 0.3212 -0.0211 -0.3094 -0.2922 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.1602B+0.9597B^2    0.6044+-0.8226i      0.9796       0.1491
## 1+0.4304B+0.8779B^2   -0.2451+-1.0387i      0.9370       0.2869
## 1-1.6909B+0.8047B^2    1.0506+-0.3727i      0.8970       0.0542
## 1+1.5891B+0.6861B^2   -1.1580+-0.3413i      0.8283       0.4544
## 1+0.9816B+0.6281B^2   -0.7814+-0.9907i      0.7925       0.3563
##   
## 

## [1] "Window ASE: 0.00056630865776623"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

Neural Network

The neural network model performs fairly well, and picks up on the weekly trend, but regresses towards a mean, and is edged out by the AR(15) model from earlier.

ts_la <- ts(la$positivity[1:239], start = '1')
x = mlp(ts_la, lags = 7, hd.auto.type = 'cv', reps=10)
plot(x)

preds <- predict(x, 7)
plot(preds)

eval_model(la$positivity,preds$mean,7,model_name = 'NNFOR', AIC_val = 0) #ASE = .000766
## Warning: attributes are not identical across measure variables;
## they will be dropped

ts_la <- ts(la$positivity[1:156], start = '1')
x = mlp(ts_la, lags = 7, m = 1, hd.auto.type = 'cv')
plot(x)

preds <- predict(x, 90)
plot(preds)

eval_model(la$positivity,preds$mean,90,model_name = 'NNFOR', AIC_val = 0) #ASE = .000852
## Warning: attributes are not identical across measure variables;
## they will be dropped

rolling <- rolling_ASE(la, e, s=7, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000637
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000982985421197367"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000125292161171138"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.000157626300527327"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000171504682305925"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.00158874760031172"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000544659687271713"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

United States

Plot the Data:

# Plot the data
plotts.sample.wge(us$positivity)

## $autplt
##  [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
##  [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  16.928178  17.041622  10.954209   3.357603   1.069112  -3.238369
##   [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
##  [13]  -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
##  [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375  -7.015352
##  [25]  -8.700809  -8.412062  -7.545175  -6.229323 -11.252124 -16.095149
##  [31] -15.372565 -27.552336 -10.346212  -9.057794 -11.461366  -7.109201
##  [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
##  [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
##  [49] -12.990975 -10.740763 -14.551855  -9.308655  -7.444421  -5.100658
##  [55]  -8.399431  -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
##  [61] -13.293199  -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
##  [67]  -8.756403 -15.569710 -11.708858 -22.389341  -7.883153 -18.081602
##  [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
##  [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
##  [85] -14.128733 -13.936455  -8.979219 -13.464535 -15.769438 -12.144105
##  [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174  -9.903788
##  [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014  -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884  -8.698403
## [121] -17.478083 -17.827845 -23.241743
## 
## $dbz
##   [1]  12.8975742  12.6104343  12.1288216  11.4480751  10.5614946   9.4601976
##   [7]   8.1330155   6.5665777   4.7459620   2.6568532   0.2916097  -2.3347667
##  [13]  -5.1453437  -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
##  [19] -11.1158574 -10.4346224  -9.7877324  -9.2573489  -8.8845067  -8.6785829
##  [25]  -8.6271903  -8.7045894  -8.8785406  -9.1161191  -9.3884024  -9.6733156
##  [31]  -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
##  [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
##  [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
##  [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
##  [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
##  [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
##  [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
##  [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
##  [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
##  [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
##  [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
##  [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607
# Look deeper at spectral density
parzen.wge(us$positivity, trun=100)
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $pzgram
##   [1]  15.382933  14.588356  12.831036   9.869648   5.787060   1.193729
##   [7]  -3.198010  -7.039404  -9.463870 -10.954884 -11.649927 -11.146329
##  [13] -10.497311 -10.387322 -10.996468 -12.241666 -14.051773 -16.357297
##  [19] -17.782568 -16.624195 -14.478599 -12.431116 -10.590963  -9.053923
##  [25]  -7.897040  -7.190547  -7.099666  -7.820660  -9.477927 -12.016628
##  [31] -14.502154 -14.318422 -11.893364  -9.673811  -8.566444  -8.883526
##  [37] -10.680694 -13.468736 -15.662532 -16.381447 -16.513312 -16.415800
##  [43] -16.292762 -16.362416 -16.490077 -16.180717 -15.264532 -14.173466
##  [49] -13.172032 -12.002961 -10.387418  -8.753405  -7.762275  -7.720243
##  [55]  -8.641456 -10.426594 -12.887273 -15.145752 -15.232845 -13.814286
##  [61] -12.874594 -12.869708 -13.340726 -13.353437 -12.555602 -11.409945
##  [67] -10.418409  -9.928099  -9.987950 -10.427637 -11.256579 -12.668644
##  [73] -14.277216 -14.778950 -14.334108 -14.338547 -14.997710 -15.512866
##  [79] -15.646402 -16.183486 -17.184966 -17.379811 -16.333944 -15.008996
##  [85] -13.685617 -12.594381 -12.257052 -12.870686 -13.936525 -14.509900
##  [91] -14.398784 -13.976794 -13.305756 -12.477274 -11.806977 -11.561958
##  [97] -11.724854 -12.038907 -12.361702 -12.843276 -13.613715 -14.408198
## [103] -14.402867 -13.168262 -11.604517 -10.672442 -10.804908 -12.107197
## [109] -14.317964 -16.503509 -17.646159 -18.062309 -18.410008 -18.446313
## [115] -17.213576 -15.075943 -13.043248 -11.492666 -10.719546 -11.075606
## [121] -12.778110 -15.616148 -17.533541

There is surpsingly little evidence of serial correlation in the US data compared to Louisiana. There are very weak peaks in the spectral density plots and a very slowly damping behavior in the autocorrelation function plots. Increasing the truncation point of the Parzen function doesn’t show significant points to indicate a significant serial correlation.

Evaluate Stationarity:

There are few indications that suggest the data is not stationary, it appears that this could be a stationary time series with some wandering behavior.

Model Construction:

Differencing the data once shows no significant cyclical pattern in the ACF.

Identify correlation structures:

Here we are going to look at different filters to try to model out some of the correlation in the data. We’ll also look at a factor table to see if we can identify a weekly or monthly trend.

# Difference the data once
us.d1 <- artrans.wge(us$positivity,1)

# Difference it again, no improvement
artrans.wge(us.d1,1)

##   [1] -0.01425418201  0.01978404489  0.00344009601 -0.06802340221  0.07352438226
##   [6]  0.00692990610 -0.00714591876 -0.02084688094 -0.00243451668  0.04964910774
##  [11] -0.03317545244 -0.05925122282  0.04131133599  0.02405539404  0.00601164353
##  [16] -0.01817377270  0.03083149410 -0.08610403374  0.08409039403 -0.02307959862
##  [21] -0.01642926846  0.01098569962 -0.09811599954  0.15692527231 -0.08792327276
##  [26]  0.03274987079 -0.01386236757  0.00918935394  0.00950420531 -0.01219184870
##  [31] -0.02247301575  0.00493196641 -0.01035396435  0.06647902814 -0.06646866380
##  [36]  0.04033066349 -0.03339489118  0.02054824224 -0.00647801193 -0.01504550299
##  [41] -0.05322120943  0.13931525485 -0.07895743252 -0.00049040986  0.01314512699
##  [46] -0.01923457995  0.02586767756 -0.02425491425  0.02450886876 -0.00465228139
##  [51] -0.01169959290  0.01236632704 -0.02144306358  0.00542332225  0.02043214764
##  [56] -0.01299864507  0.00739442008 -0.00979988508  0.00423135860 -0.02960535336
##  [61]  0.05470381881 -0.02719446153  0.00982324167 -0.01199966873  0.00805208606
##  [66] -0.00623128961  0.00376016574 -0.00321244605  0.00046827232  0.00967439924
##  [71] -0.00039847540 -0.01103591006  0.00557147233 -0.00879606612  0.01939398534
##  [76] -0.01744094841  0.00899981221 -0.00453842682  0.00960374829 -0.00902539121
##  [81]  0.00071623524 -0.00215174346  0.00057748709  0.00145547567 -0.00053317682
##  [86]  0.00332735306  0.00265811119 -0.00623099656 -0.00110895510  0.01109982666
##  [91] -0.00786465370 -0.00713998322  0.01679738401 -0.01340796604  0.00140554934
##  [96]  0.00886218895 -0.00864097890  0.00711356390 -0.00487214075  0.00488677486
## [101] -0.00584010118  0.00130517473 -0.00229574827  0.02632837457 -0.04012080831
## [106]  0.01463289507  0.00921421538 -0.00786276163 -0.00956919410  0.02247511852
## [111] -0.00867981153 -0.00477364349 -0.00200677029  0.01103878341 -0.02503295136
## [116]  0.01482554512  0.01315600697  0.00117340963 -0.02472644128  0.01051180247
## [121]  0.00847784399 -0.01529062714  0.01004135529 -0.00054810898 -0.00068255160
## [126] -0.00003667196  0.00499935746 -0.01283051447  0.00653668209  0.00001297137
## [131]  0.00435733193 -0.00366747653  0.00009589148 -0.00386304015  0.00230426213
## [136] -0.00426560014  0.00520113038  0.00308756001 -0.00116121005  0.00533521705
## [141] -0.01614762899  0.01470606243 -0.01621035783  0.01068486672  0.00547535911
## [146] -0.00548022191 -0.00180068959  0.00571874086 -0.00010420098 -0.01070456362
## [151] -0.00465714413  0.02870969236  0.01514774804 -0.08250772757  0.05829342260
## [156] -0.00198719651 -0.02080531661  0.01639823857  0.00370790673 -0.00407842935
## [161] -0.00434054222  0.00265285183 -0.00085254569  0.00250766698 -0.00625128498
## [166]  0.00942351027  0.00338309003 -0.01553426756  0.00769235953 -0.00118990894
## [171] -0.00037389962 -0.00315165280  0.01446704665 -0.02152086746  0.02277328413
## [176] -0.01200183966 -0.00177860185 -0.00346617760  0.00911724230 -0.00659330821
## [181]  0.01005544004 -0.00283078439 -0.00248659228 -0.00708970815  0.00880266567
## [186]  0.02081367549 -0.04811764792  0.03040237498 -0.00930103949  0.00235031086
## [191] -0.00360694153  0.00316510150  0.01613333560 -0.00720811766 -0.02667476139
## [196]  0.02312455662  0.00430214086 -0.01507066638  0.00371924998 -0.00424231109
## [201]  0.01639379996  0.00058917864 -0.01267489139 -0.00217868997  0.01234496297
## [206] -0.01613330064  0.01052661225  0.00344839781  0.00460398622 -0.01318871581
## [211]  0.00352505023 -0.00084230457 -0.00159297749  0.00027587920  0.00862199954
## [216]  0.00155555540 -0.00407252270 -0.00513442594 -0.00960502040  0.01307353256
## [221]  0.00264711510  0.00880620228 -0.01204944398 -0.00666575615  0.00578262594
## [226]  0.00311978586 -0.01173173304  0.00533528809  0.01865195916 -0.01111606447
## [231] -0.00982475902  0.01001397542 -0.01088472895  0.00645926430 -0.00394642266
## [236]  0.01856996819 -0.00723121433 -0.01788513476  0.01628207270  0.00493381307
## [241] -0.01441135849 -0.00594183904  0.03534821567 -0.02865985992
# Try various seasonal models
us.s7 <- artrans.wge(us$positivity, c(rep(0,6),1))

us.s12 <- artrans.wge(us$positivity, c(rep(0,11),1))

us.d1.s7 <- artrans.wge(us.d1, c(rep(0,6),1))

us.d1.s12 <- artrans.wge(us.d1, c(rep(0,11),1))

# Overfit table
est.ar.wge(us$positivity, p=15)
## 
## Coefficients of Original polynomial:  
## 0.3866 0.1849 0.1179 0.1508 0.1336 -0.0545 0.1364 0.2048 0.0220 0.0994 -0.1207 -0.1809 -0.0703 0.1205 -0.1875 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9642B+0.9671B^2    1.0155+-0.0524i      0.9834       0.0082
## 1-0.3434B+0.8525B^2    0.2014+-1.0642i      0.9233       0.2202
## 1+0.9180B             -1.0893               0.9180       0.5000
## 1+0.4309B+0.8427B^2   -0.2557+-1.0589i      0.9180       0.2877
## 1-1.3843B+0.8248B^2    0.8392+-0.7129i      0.9082       0.1121
## 1+1.2323B+0.8146B^2   -0.7564+-0.8096i      0.9025       0.3696
## 1+1.6110B+0.8080B^2   -0.9970+-0.4937i      0.8989       0.4268
## 1-0.8868B+0.5415B^2    0.8189+-1.0845i      0.7358       0.1471
##   
##   
## 
## Coefficients of Original polynomial:  
## 0.3866 0.1849 0.1179 0.1508 0.1336 -0.0545 0.1364 0.2048 0.0220 0.0994 -0.1207 -0.1809 -0.0703 0.1205 -0.1875 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9642B+0.9671B^2    1.0155+-0.0524i      0.9834       0.0082
## 1-0.3434B+0.8525B^2    0.2014+-1.0642i      0.9233       0.2202
## 1+0.9180B             -1.0893               0.9180       0.5000
## 1+0.4309B+0.8427B^2   -0.2557+-1.0589i      0.9180       0.2877
## 1-1.3843B+0.8248B^2    0.8392+-0.7129i      0.9082       0.1121
## 1+1.2323B+0.8146B^2   -0.7564+-0.8096i      0.9025       0.3696
## 1+1.6110B+0.8080B^2   -0.9970+-0.4937i      0.8989       0.4268
## 1-0.8868B+0.5415B^2    0.8189+-1.0845i      0.7358       0.1471
##   
## 
## $phi
##  [1]  0.38657693  0.18487981  0.11787070  0.15076990  0.13357768 -0.05449747
##  [7]  0.13638095  0.20481227  0.02200580  0.09936180 -0.12065274 -0.18089166
## [13] -0.07032510  0.12048686 -0.18746678
## 
## $res
##   [1]  0.01658565790  0.01375283230 -0.00090975462  0.00845540202  0.01787093204
##   [6] -0.03984666673 -0.00642488247  0.01431184179  0.01779142314  0.01080242259
##  [11]  0.00663355927  0.04102041625  0.03926094547 -0.01183329929 -0.01332829485
##  [16]  0.01014091369  0.01845543645  0.00884849996  0.04057279959 -0.02456491855
##  [21]  0.01134078334  0.01787148894  0.00870671436  0.02431789715 -0.07283238027
##  [26] -0.00517819506 -0.03437950208  0.00118715503 -0.01032289425  0.00644319619
##  [31]  0.00652596746  0.01987781181  0.00488724220  0.00373016664 -0.01721744312
##  [36]  0.02418732184 -0.01049618328  0.01171728082  0.00293030831 -0.01518432094
##  [41] -0.00155289996 -0.01039100053 -0.06795034138  0.02871827380  0.00452748430
##  [46] -0.01239578134  0.00854411475 -0.00827701124  0.00216160779 -0.00435294323
##  [51]  0.02690416243  0.01475487501  0.01892083537  0.00496840193 -0.01299866922
##  [56] -0.00935331316  0.01502989285 -0.00836566651  0.00641615521  0.00149736569
##  [61] -0.00310893443 -0.03220818092  0.00693531442  0.00508118048  0.00916087672
##  [66]  0.00403499757  0.00518738405 -0.00595027914  0.00420244633  0.00369506342
##  [71] -0.00360685454  0.00856749870  0.00556391548 -0.00256806545 -0.00069150154
##  [76] -0.00312190783  0.00158672696 -0.00127510904  0.00152979208 -0.00120194631
##  [81]  0.00433291352 -0.00098428801 -0.00156131332 -0.00254140231 -0.00486746293
##  [86] -0.00545535807 -0.00710667225 -0.00227452103  0.00061060139 -0.00133946919
##  [91] -0.00769660571  0.00591309232  0.00306698063 -0.00369564838  0.00718262066
##  [96] -0.00001671539 -0.00461294328  0.00401196919 -0.00159613819  0.00109499152
## [101]  0.00196209263  0.00371501321 -0.00167261727  0.00186236438 -0.00299427170
## [106]  0.01941514673 -0.00382930045 -0.00384991841  0.00190630333 -0.00071041088
## [111] -0.01224783761  0.00632026187  0.00502788605 -0.00259796971 -0.00119884956
## [116]  0.00344108999 -0.01517812106 -0.00822784957  0.00497434318  0.00772444255
## [121] -0.00229308135 -0.00430162836 -0.00085720848 -0.00944426026  0.00157793679
## [126]  0.00018534116 -0.00056310649 -0.00126815864  0.00390301921 -0.00977286373
## [131] -0.00092462071 -0.00053444002 -0.00068972220 -0.00069105253  0.00380122262
## [136] -0.00448148251 -0.00643582696 -0.00728575132 -0.00634107135 -0.00040158021
## [141]  0.00118932466  0.00684192504 -0.00612621632  0.00392189653 -0.00818374001
## [146] -0.00325862021  0.00308286413  0.00267336077 -0.00223351023  0.00232775916
## [151]  0.00372790022 -0.00494997796 -0.01113158303  0.00976643002  0.03693895275
## [156] -0.02817729817 -0.00566048758 -0.00240558087 -0.01928741546 -0.00913914196
## [161]  0.00412288863 -0.00488766390 -0.00643902903 -0.00091385188 -0.00827121966
## [166]  0.00636653750  0.00466488156  0.00003018588  0.00509632278  0.00592229541
## [171] -0.00489662044 -0.00466327616 -0.00237286795 -0.00815518027  0.00495705805
## [176] -0.01135691942  0.00368596891  0.00028501350 -0.00433408475 -0.00817656855
## [181]  0.00080322115 -0.00704205020  0.00114529007  0.00694825730  0.00090665491
## [186] -0.00526555278 -0.00078038882  0.02232359024 -0.01211795724  0.00344900910
## [191] -0.00690745033 -0.00855145879 -0.01250710823 -0.00693007971  0.00369252631
## [196]  0.00846396477 -0.01386141770 -0.00732336885  0.00637722556 -0.00252184851
## [201] -0.00431535013 -0.01269272632  0.00126723235  0.00387963806  0.00032530003
## [206] -0.00700880871  0.00587812304 -0.00701177611 -0.00660973973  0.00352045916
## [211]  0.01335070103 -0.00109421533 -0.00276350971 -0.00446999423 -0.00667799040
## [216] -0.00394059592  0.00040448468  0.00528969534  0.00592415760  0.00070268984
## [221] -0.01447748759 -0.00501264793  0.00045067434  0.01262251779  0.00709674750
## [226] -0.00096957046 -0.00472347849 -0.00187216773 -0.00895473279 -0.00658333715
## [231]  0.01023025936  0.00504505014 -0.00489649435  0.00159401036 -0.00531932415
## [236] -0.00433809527 -0.00418639105  0.00973786582  0.01178895197 -0.00241604198
## [241]  0.00069451250  0.00853576955  0.00406154919 -0.00520709593  0.02140959751
## [246]  0.00587643257
## 
## $avar
## [1] 0.0001537157
## 
## $aic
## [1] -8.650324
## 
## $aicc
## [1] -7.631283
## 
## $bic
## [1] -8.422335
# Generate a factor table for a (1-B^7)
tswge::factor.wge(c(rep(0,6),1))
## 
## Coefficients of Original polynomial:  
## 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.0000B              1.0000               1.0000       0.0000
## 1+0.4450B+1.0000B^2   -0.2225+-0.9749i      1.0000       0.2857
## 1-1.2470B+1.0000B^2    0.6235+-0.7818i      1.0000       0.1429
## 1+1.8019B+1.0000B^2   -0.9010+-0.4339i      1.0000       0.4286
##   
## 

Estimate parameters and check residuals:

The ACF and PACF plots look like white noise, as expected. The Ljung-Box test fails to reject the null hyphothesis that the residuals are white noise. The generated spectral densities look like they have a fairly good fit, the final peak is overpronounced in our generations compared to the actual data. The ACF comparisons are all over the board, but with peaks in the same general spots (lags at 7, 14, 21). The generated realizations look quite a bit different, but the spectral densities and ACF plots look fairly similar, indicating this could be a useful model.

# Estimate the model params
aic5.wge(us$positivity, p=0:20, q=0:2, type='bic') # aic = 19,1 > 18,0 : bic = 1,1 > 9,1 > 2,1
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Error in aic calculation at 20 1 
## Five Smallest Values of  bic
##       p    q        bic
## 5     1    1  -8.487898
## 29    9    1  -8.468857
## 8     2    1  -8.465559
## 6     1    2  -8.465077
## 35   11    1  -8.462634
e <- est.arma.wge(us$positivity, p=9, q=1)
## 
## Coefficients of Original polynomial:  
## 1.2711 -0.1631 -0.0699 0.0662 0.0028 -0.1636 0.1612 0.1144 -0.2295 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9861B+0.9886B^2    1.0045+-0.0501i      0.9943       0.0079
## 1-0.0326B+0.7214B^2    0.0226+-1.1772i      0.8493       0.2469
## 1-1.1925B+0.6943B^2    0.8588+-0.8384i      0.8332       0.1231
## 1+1.2068B+0.6319B^2   -0.9548+-0.8190i      0.7949       0.3872
## 1+0.7333B             -1.3636               0.7333       0.5000
##   
## 
# Ljung Box Test still shows white noise residuals
ljung.wge(e$res)
## Obs -0.03356559 -0.04604162 -0.004168052 0.01547628 -0.01688291 -0.03271148 -0.0009237203 0.003706762 0.06424109 0.1178936 -0.03902071 -0.1007295 -0.06082771 0.1679569 -0.04072279 -0.1047464 -0.02937219 0.1135898 0.1318257 0.02116977 0.05601096 -0.05664106 0.07570874 -0.1847351
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 24
## 
## $chi.square
## [1] 41.79917
## 
## $df
## [1] 24
## 
## $pval
## [1] 0.0135891
ljung.wge(e$res, K = 48)
## Obs -0.03356559 -0.04604162 -0.004168052 0.01547628 -0.01688291 -0.03271148 -0.0009237203 0.003706762 0.06424109 0.1178936 -0.03902071 -0.1007295 -0.06082771 0.1679569 -0.04072279 -0.1047464 -0.02937219 0.1135898 0.1318257 0.02116977 0.05601096 -0.05664106 0.07570874 -0.1847351 -0.04456868 -0.1097331 -0.008943643 0.1029134 0.03469855 -0.1181843 -0.04739403 0.0068443 0.01577266 -0.04474552 0.08704028 0.02497461 0.1092873 -0.08576074 0.01382154 -0.01940963 0.0472867 0.06192417 -0.01025361 -0.02817321 -0.0161036 0.04494596 -0.02163542 0.009431244
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 48
## 
## $chi.square
## [1] 65.38729
## 
## $df
## [1] 48
## 
## $pval
## [1] 0.04816801
acf(e$res)

pacf(e$res)

dev.off()
## null device 
##           1
#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(us$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)

for( i in 1: sims)
{
   SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 0, phi = e$phi, plot ="FALSE"), plot = "FALSE")
   lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}

#Compare ACFs
sims = 5
ACF = acf(us$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)

for( i in 1: sims)
{
   ACF2 = acf(gen.aruma.wge(246, s = 0, phi = e$phi, plot = "FALSE"), plot = "FALSE")
   lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}

#Compare Generated Realizations 
eGen = gen.aruma.wge(246, s = 0, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)
## $autplt
##  [1] 1.0000000 0.9935339 0.9838052 0.9713461 0.9562024 0.9380209 0.9171481
##  [8] 0.8939896 0.8687276 0.8414615 0.8122268 0.7814705 0.7491547 0.7152164
## [15] 0.6803175 0.6438795 0.6061015 0.5667866 0.5267571 0.4858514 0.4440230
## [22] 0.4014926 0.3584973 0.3151003 0.2707877 0.2260548
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]   4.785980  20.460683   6.878691  -9.792577  -2.442456  -6.769525
##   [7]  -2.114508  -5.513326  -4.666739  -8.681357  -6.524804  -7.989461
##  [13]  -9.545319 -10.547635  -9.767938 -13.914741 -13.222626 -11.588570
##  [19] -14.143523 -12.485914 -14.992938 -14.761362 -19.147332 -12.042511
##  [25] -18.455932 -16.026129 -15.447145 -20.189811 -18.951746 -35.972499
##  [31] -15.036756 -21.535024 -14.861839 -22.526411 -18.001167 -19.586483
##  [37] -18.745639 -18.769383 -22.678378 -21.667724 -19.555406 -18.744277
##  [43] -19.948553 -19.908577 -21.884591 -23.542054 -21.372526 -20.604435
##  [49] -18.870470 -21.058358 -18.361203 -18.387750 -22.869790 -22.423839
##  [55] -21.306296 -22.262206 -24.710437 -24.561472 -19.060970 -26.298557
##  [61] -22.948965 -22.721761 -21.603496 -20.555001 -22.313456 -21.712554
##  [67] -26.106076 -25.014777 -23.047509 -23.253366 -24.970398 -28.859453
##  [73] -24.665586 -23.240980 -31.541241 -23.468987 -25.489892 -24.517900
##  [79] -25.702613 -28.213259 -25.260366 -28.133953 -25.136608 -23.862811
##  [85] -22.097075 -27.666242 -24.031868 -23.117571 -27.022319 -27.565200
##  [91] -27.458709 -23.068928 -27.052203 -27.098103 -31.964080 -26.061005
##  [97] -28.751894 -26.069454 -26.578042 -27.960504 -29.763770 -26.343833
## [103] -26.969098 -23.194316 -26.482595 -25.923156 -25.220976 -31.365838
## [109] -26.645385 -23.904860 -26.940378 -28.862589 -28.254769 -27.305370
## [115] -30.382414 -24.531366 -28.769833 -23.600321 -28.620172 -26.270184
## [121] -26.089719 -23.583943 -26.787258
## 
## $dbz
##   [1]  12.925445  12.659061  12.212530  11.581967  10.761845   9.744968
##   [7]   8.522543   7.084561   5.420919   3.524292   1.397019  -0.933425
##  [13]  -3.383545  -5.769855  -7.805983  -9.251522 -10.119925 -10.619053
##  [19] -10.932032 -11.147277 -11.308923 -11.465576 -11.679084 -12.009401
##  [25] -12.498549 -13.161416 -13.979962 -14.897436 -15.816880 -16.618552
##  [31] -17.207066 -17.561524 -17.736164 -17.811584 -17.850539 -17.886983
##  [37] -17.936897 -18.012207 -18.127776 -18.299997 -18.540416 -18.848757
##  [43] -19.208424 -19.586555 -19.940482 -20.230689 -20.435419 -20.557951
##  [49] -20.621475 -20.656312 -20.688648 -20.735738 -20.806683 -20.905491
##  [55] -21.033596 -21.190547 -21.373111 -21.573992 -21.781707 -21.982598
##  [61] -22.164812 -22.322612 -22.458625 -22.582582 -22.707176 -22.843271
##  [67] -22.996590 -23.166879 -23.349321 -23.537228 -23.724637 -23.907666
##  [73] -24.084143 -24.252047 -24.407932 -24.546465 -24.661534 -24.748438
##  [79] -24.806027 -24.837644 -24.850502 -24.853978 -24.857791 -24.870730
##  [85] -24.900084 -24.951444 -25.028498 -25.132602 -25.262198 -25.412357
##  [91] -25.574816 -25.738750 -25.892239 -26.024081 -26.125357 -26.190261
##  [97] -26.216163 -26.203322 -26.154771 -26.076539 -25.977872 -25.870918
## [103] -25.769584 -25.687730 -25.637172 -25.625990 -25.657378 -25.729124
## [109] -25.833696 -25.958940 -26.089447 -26.208624 -26.301300 -26.356388
## [115] -26.368892 -26.340592 -26.279221 -26.196519 -26.105888 -26.020309
## [121] -25.950849 -25.905790 -25.890203
plotts.sample.wge(us$positivity)
## $autplt
##  [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
##  [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  16.928178  17.041622  10.954209   3.357603   1.069112  -3.238369
##   [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
##  [13]  -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
##  [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375  -7.015352
##  [25]  -8.700809  -8.412062  -7.545175  -6.229323 -11.252124 -16.095149
##  [31] -15.372565 -27.552336 -10.346212  -9.057794 -11.461366  -7.109201
##  [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
##  [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
##  [49] -12.990975 -10.740763 -14.551855  -9.308655  -7.444421  -5.100658
##  [55]  -8.399431  -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
##  [61] -13.293199  -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
##  [67]  -8.756403 -15.569710 -11.708858 -22.389341  -7.883153 -18.081602
##  [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
##  [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
##  [85] -14.128733 -13.936455  -8.979219 -13.464535 -15.769438 -12.144105
##  [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174  -9.903788
##  [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014  -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884  -8.698403
## [121] -17.478083 -17.827845 -23.241743
## 
## $dbz
##   [1]  12.8975742  12.6104343  12.1288216  11.4480751  10.5614946   9.4601976
##   [7]   8.1330155   6.5665777   4.7459620   2.6568532   0.2916097  -2.3347667
##  [13]  -5.1453437  -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
##  [19] -11.1158574 -10.4346224  -9.7877324  -9.2573489  -8.8845067  -8.6785829
##  [25]  -8.6271903  -8.7045894  -8.8785406  -9.1161191  -9.3884024  -9.6733156
##  [31]  -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
##  [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
##  [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
##  [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
##  [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
##  [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
##  [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
##  [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
##  [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
##  [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
##  [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
##  [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607
# Check performance
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(9,1)',AIC_val = e$aic) #ASE .000121

preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(9,1)', AIC_val = e$aic) #ASE .00115

rolling <- rolling_ASE(us, e, s=0, horizon=7, model_name = 'AR(9,1)' ) #ASE .000096
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.00019584690928318"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.0000286092864906074"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000184254204829249"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000264139426854921"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.0000280660276168899"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.0000435508295999369"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).

Fitting additional models:

ARMA(18,1) with a weekly trend

The ARMA(18,1) with a weekly trend shows a slight drop in performance in the 7 day prediction, but much better performance in the 90 day trend. It is the most balanced model, having the best performance in 90 day predictions and competitive performance in the rolling window ASE.

aic5.wge(us.s7, p=0:20, q=0:5, type='bic') #aic = 20,1 > 20,2 > 18,1 : bic = 20,1 > 18,1
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Error in aic calculation at 13 1 
## Error in aic calculation at 15 1 
## Error in aic calculation at 15 4 
## Error in aic calculation at 15 5 
## Error in aic calculation at 16 3 
## Five Smallest Values of  bic
##        p    q        bic
## 122   20    1  -8.432983
## 110   18    1  -8.421151
## 116   19    1  -8.403711
## 123   20    2  -8.397647
## 53     8    4  -8.354801
e <- est.arma.wge(us.s7, p=18, q=1)
## 
## Coefficients of Original polynomial:  
## 1.2038 -0.1377 0.0970 -0.1074 -0.0016 -0.0903 -0.6689 1.0683 -0.2391 0.2188 -0.4459 -0.0299 0.0073 -0.1330 0.4938 -0.2741 0.3459 -0.3470 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9871B+0.9916B^2    1.0019+-0.0674i      0.9958       0.0107
## 1-1.5531B+0.9722B^2    0.7988+-0.6249i      0.9860       0.1057
## 1+1.9252B+0.9394B^2   -1.0247+-0.1205i      0.9692       0.4814
## 1+1.1191B+0.9092B^2   -0.6154+-0.8492i      0.9535       0.3498
## 1+1.4737B+0.9061B^2   -0.8132+-0.6651i      0.9519       0.3909
## 1-1.7816B+0.9039B^2    0.9855+-0.3675i      0.9507       0.0568
## 1-0.3277B+0.8822B^2    0.1857+-1.0484i      0.9392       0.2221
## 1-0.5258B+0.8490B^2    0.3096+-1.0402i      0.9214       0.2040
## 1+0.4535B+0.6870B^2   -0.3301+-1.1604i      0.8289       0.2941
##   
## 
# Ljung Box Test still shows white noise residuals
ljung.wge(artrans.wge(us.s7, phi.tr = e$phi))
## Obs -0.5065971 -0.003026183 0.002316859 0.08036559 -0.1209084 0.07023771 0.04574104 -0.03331472 0.03382918 -0.103711 0.001563584 0.07111938 0.1553696 -0.2417995 -0.04760502 0.2130846 -0.1512536 0.01712802 0.1140875 0.00473265 -0.1590841 0.1094058 -0.09090864 0.1008303
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 24
## 
## $chi.square
## [1] 121.815
## 
## $df
## [1] 24
## 
## $pval
## [1] 0.000000000000004551914
ljung.wge(artrans.wge(us.s7, phi.tr = e$phi), K = 48)

## Obs -0.5065971 -0.003026183 0.002316859 0.08036559 -0.1209084 0.07023771 0.04574104 -0.03331472 0.03382918 -0.103711 0.001563584 0.07111938 0.1553696 -0.2417995 -0.04760502 0.2130846 -0.1512536 0.01712802 0.1140875 0.00473265 -0.1590841 0.1094058 -0.09090864 0.1008303 0.02028824 -0.03678615 -0.05092357 0.02339115 0.0240892 0.01927424 -0.04905638 0.04213179 -0.0005581815 -0.1033197 0.1318134 -0.09368698 0.08528423 -0.02169268 -0.03422522 -0.04499843 0.08258557 -0.03179607 -0.01507485 0.02660707 -0.02762987 0.0711119 -0.08979131 -0.008203965
## $test
## [1] "Ljung-Box test"
## 
## $K
## [1] 48
## 
## $chi.square
## [1] 143.4359
## 
## $df
## [1] 48
## 
## $pval
## [1] 0.00000000001924017
acf(e$res)
pacf(e$res)
dev.off()
## null device 
##           1
#Compare Spectral Densities
sims = 5
SpecDen = parzen.wge(us$positivity, plot = "FALSE")
plot(SpecDen$freq, SpecDen$pzgram, type = "l", lwd = 6)

for( i in 1: sims)
{
   SpecDen2 = parzen.wge(gen.aruma.wge(246,s = 7, phi = e$phi, plot ="FALSE"), plot = "FALSE")
   lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}

#Compare ACFs
sims = 5
ACF = acf(us$positivity, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)

for( i in 1: sims)
{
   ACF2 = acf(gen.aruma.wge(246, s = 7, phi = e$phi, plot = "FALSE"), plot = "FALSE")
   lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}

#Compare Generated Realizations 
eGen = gen.aruma.wge(246, s = 7, phi = e$phi, vara = e$avar)
plotts.sample.wge(eGen)
## $autplt
##  [1] 1.0000000 0.9848470 0.9714933 0.9750294 0.9660213 0.9445656 0.9409767
##  [8] 0.9406086 0.9164590 0.8947200 0.8887335 0.8694487 0.8382897 0.8244560
## [15] 0.8136619 0.7805097 0.7512321 0.7387271 0.7140631 0.6782452 0.6590084
## [22] 0.6421241 0.6032096 0.5685146 0.5501356 0.5211394
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  17.4883412  17.8759163   3.3192846  -3.8525491  -0.3197144 -21.6714749
##   [7] -22.8063713 -10.2392453 -12.0673705 -12.2214139 -13.9805703 -18.7730531
##  [13] -18.8355978 -14.2310658 -18.3550959 -24.3257081 -18.1713112 -18.9508838
##  [19] -21.0905082 -22.5490074 -22.9798514 -27.1738870 -22.4637158 -26.0337988
##  [25] -16.3206084  -8.5546634 -15.3803215 -19.2883534 -17.8332777 -18.1899535
##  [31] -26.4043277 -23.8962853 -22.7948280 -16.8336539  -6.1759350 -17.3886878
##  [37] -26.3712238 -30.6911752 -30.0931877 -25.9624947 -23.5436142 -25.9359153
##  [43] -25.3700146 -23.9993226 -25.8558264 -26.4139571 -23.7900585 -30.4384324
##  [49] -30.0297904 -27.3820589 -30.1075893 -19.9753587 -21.7443303 -19.2880151
##  [55] -26.0065533 -34.3012107 -21.6232692 -24.7615997 -27.7416393 -26.4778644
##  [61] -26.6516604 -21.3567198 -22.4913621 -21.6969354 -20.1932738 -20.4355328
##  [67] -17.8545927 -16.6305613 -12.4553796  -2.0051496  -6.0649279 -18.2027733
##  [73] -22.1637145 -21.1998955 -24.6167114 -27.6646930 -30.7627963 -28.9387210
##  [79] -37.0841490 -35.0228031 -33.7229017 -30.6743821 -35.4618060 -32.0714755
##  [85] -30.0881307 -29.4898690 -37.9232247 -40.4322393 -36.4403348 -37.8417324
##  [91] -31.3156666 -36.6701417 -32.4208118 -37.4804681 -31.6757217 -28.1692244
##  [97] -30.6926889 -37.2504701 -43.4538675 -29.2992604 -34.3024720 -31.7441604
## [103] -28.5361267 -22.1853685 -24.1843615 -18.4839851 -27.5289666 -33.4879999
## [109] -39.8646199 -33.4828829 -34.6946852 -38.9555566 -35.0056370 -46.1802585
## [115] -38.4655865 -34.6454560 -39.9518903 -36.9388807 -24.9050466 -41.8672693
## [121] -39.9673409 -38.8068328 -34.7290668
## 
## $dbz
##   [1]  13.1505985  12.8508351  12.3476953  11.6357495  10.7071467   9.5513936
##   [7]   8.1551522   6.5022292   4.5742086   2.3529001  -0.1723734  -2.9851629
##  [13]  -5.9943454  -8.9461610 -11.3785406 -12.8764690 -13.4514878 -13.3828185
##  [19] -12.9305630 -12.3235211 -11.7616150 -11.3810556 -11.2483531 -11.3764193
##  [25] -11.7386223 -12.2740811 -12.8899611 -13.4736343 -13.9256735 -14.2014856
##  [31] -14.3245535 -14.3576413 -14.3642960 -14.3912229 -14.4710836 -14.6318058
##  [37] -14.9027654 -15.3154018 -15.8998317 -16.6793579 -17.6628684 -18.8327983
##  [43] -20.1258665 -21.4114139 -22.4978380 -23.2156270 -23.5363512 -23.5726581
##  [49] -23.4591612 -23.2742952 -23.0485995 -22.8036367 -22.5760745 -22.4147948
##  [55] -22.3601135 -22.4148600 -22.5064242 -22.4471161 -21.9635668 -20.8921253
##  [61] -19.3571507 -17.6410485 -15.9700130 -14.4625165 -13.1665084 -12.0964210
##  [67] -11.2533126 -10.6339936 -10.2348609 -10.0534452 -10.0890256 -10.3428890
##  [73] -10.8184658 -11.5214051 -12.4595315 -13.6424492 -15.0801763 -16.7792701
##  [79] -18.7327507 -20.8959192 -23.1374554 -25.1828816 -26.6789790 -27.4915621
##  [85] -27.8444184 -28.0356068 -28.2163791 -28.4034520 -28.5572936 -28.6450395
##  [91] -28.6696151 -28.6616936 -28.6512341 -28.6425116 -28.6051411 -28.4842864
##  [97] -28.2288276 -27.8241888 -27.3050207 -26.7379019 -26.1926286 -25.7235865
## [103] -25.3653040 -25.1357367 -25.0412690 -25.0807136 -25.2475968 -25.5307047
## [109] -25.9130085 -26.3693131 -26.8636964 -27.3492391 -27.7738304 -28.0941293
## [115] -28.2925476 -28.3849520 -28.4114379 -28.4166945 -28.4327401 -28.4708364
## [121] -28.5224704 -28.5672362 -28.5849455
plotts.sample.wge(us$positivity)
## $autplt
##  [1] 1.0000000 0.9470164 0.9341278 0.9260111 0.9198625 0.9088272 0.8909389
##  [8] 0.8836047 0.8738064 0.8521265 0.8325090 0.8001230 0.7724089 0.7494712
## [15] 0.7330396 0.6973509 0.6643364 0.6399202 0.6205127 0.5900024 0.5527729
## [22] 0.5218086 0.4844797 0.4575403 0.4134570 0.3899807
## 
## $freq
##   [1] 0.004065041 0.008130081 0.012195122 0.016260163 0.020325203 0.024390244
##   [7] 0.028455285 0.032520325 0.036585366 0.040650407 0.044715447 0.048780488
##  [13] 0.052845528 0.056910569 0.060975610 0.065040650 0.069105691 0.073170732
##  [19] 0.077235772 0.081300813 0.085365854 0.089430894 0.093495935 0.097560976
##  [25] 0.101626016 0.105691057 0.109756098 0.113821138 0.117886179 0.121951220
##  [31] 0.126016260 0.130081301 0.134146341 0.138211382 0.142276423 0.146341463
##  [37] 0.150406504 0.154471545 0.158536585 0.162601626 0.166666667 0.170731707
##  [43] 0.174796748 0.178861789 0.182926829 0.186991870 0.191056911 0.195121951
##  [49] 0.199186992 0.203252033 0.207317073 0.211382114 0.215447154 0.219512195
##  [55] 0.223577236 0.227642276 0.231707317 0.235772358 0.239837398 0.243902439
##  [61] 0.247967480 0.252032520 0.256097561 0.260162602 0.264227642 0.268292683
##  [67] 0.272357724 0.276422764 0.280487805 0.284552846 0.288617886 0.292682927
##  [73] 0.296747967 0.300813008 0.304878049 0.308943089 0.313008130 0.317073171
##  [79] 0.321138211 0.325203252 0.329268293 0.333333333 0.337398374 0.341463415
##  [85] 0.345528455 0.349593496 0.353658537 0.357723577 0.361788618 0.365853659
##  [91] 0.369918699 0.373983740 0.378048780 0.382113821 0.386178862 0.390243902
##  [97] 0.394308943 0.398373984 0.402439024 0.406504065 0.410569106 0.414634146
## [103] 0.418699187 0.422764228 0.426829268 0.430894309 0.434959350 0.439024390
## [109] 0.443089431 0.447154472 0.451219512 0.455284553 0.459349593 0.463414634
## [115] 0.467479675 0.471544715 0.475609756 0.479674797 0.483739837 0.487804878
## [121] 0.491869919 0.495934959 0.500000000
## 
## $db
##   [1]  16.928178  17.041622  10.954209   3.357603   1.069112  -3.238369
##   [7] -11.427364 -14.831916 -13.181967 -26.879328 -10.640049 -13.104348
##  [13]  -6.897898 -10.769553 -12.587656 -16.736482 -15.880304 -19.355323
##  [19] -19.864320 -18.877313 -12.131779 -10.974583 -16.180375  -7.015352
##  [25]  -8.700809  -8.412062  -7.545175  -6.229323 -11.252124 -16.095149
##  [31] -15.372565 -27.552336 -10.346212  -9.057794 -11.461366  -7.109201
##  [37] -14.240281 -13.075078 -23.877820 -15.330952 -14.475847 -20.765059
##  [43] -13.129351 -25.189770 -19.078147 -17.562843 -19.295449 -16.258703
##  [49] -12.990975 -10.740763 -14.551855  -9.308655  -7.444421  -5.100658
##  [55]  -8.399431  -9.881040 -17.975943 -17.462588 -17.722631 -13.792923
##  [61] -13.293199  -9.316456 -15.869391 -14.012814 -14.465016 -12.635251
##  [67]  -8.756403 -15.569710 -11.708858 -22.389341  -7.883153 -18.081602
##  [73] -16.077320 -13.828997 -15.999068 -10.429352 -26.201640 -19.135996
##  [79] -14.542536 -16.213531 -16.985225 -20.728879 -15.016682 -31.829738
##  [85] -14.128733 -13.936455  -8.979219 -13.464535 -15.769438 -12.144105
##  [91] -14.497446 -17.272748 -12.667286 -11.379894 -13.465174  -9.903788
##  [97] -13.261966 -12.491846 -11.986603 -11.921460 -13.694377 -12.698553
## [103] -16.097000 -15.677901 -11.088504 -11.527014  -9.506795 -12.681243
## [109] -14.748674 -18.413266 -16.493124 -18.659561 -18.956181 -17.539240
## [115] -23.088025 -19.681462 -12.465460 -11.043980 -10.130884  -8.698403
## [121] -17.478083 -17.827845 -23.241743
## 
## $dbz
##   [1]  12.8975742  12.6104343  12.1288216  11.4480751  10.5614946   9.4601976
##   [7]   8.1330155   6.5665777   4.7459620   2.6568532   0.2916097  -2.3347667
##  [13]  -5.1453437  -7.9215132 -10.2286515 -11.5956090 -11.9705143 -11.6983430
##  [19] -11.1158574 -10.4346224  -9.7877324  -9.2573489  -8.8845067  -8.6785829
##  [25]  -8.6271903  -8.7045894  -8.8785406  -9.1161191  -9.3884024  -9.6733156
##  [31]  -9.9562335 -10.2290872 -10.4895332 -10.7412048 -10.9945308 -11.2664674
##  [37] -11.5777577 -11.9475163 -12.3857167 -12.8839971 -13.4053244 -13.8761561
##  [43] -14.1925522 -14.2554748 -14.0264019 -13.5536343 -12.9396293 -12.2881747
##  [49] -11.6760017 -11.1508180 -10.7394607 -10.4558243 -10.3057268 -10.2889439
##  [55] -10.3991649 -10.6224555 -10.9347537 -11.2993911 -11.6668145 -11.9799663
##  [61] -12.1877268 -12.2626791 -12.2120173 -12.0727873 -11.8955891 -11.7286297
##  [67] -11.6092436 -11.5622931 -11.6019261 -11.7339670 -11.9576927 -12.2666738
##  [73] -12.6487555 -13.0854284 -13.5510287 -14.0125543 -14.4313384 -14.7679001
##  [79] -14.9901220 -15.0821781 -15.0493866 -14.9156322 -14.7148182 -14.4811756
##  [85] -14.2423645 -14.0164221 -13.8117305 -13.6288608 -13.4634734 -13.3096700
##  [91] -13.1631127 -13.0231056 -12.8930232 -12.7790349 -12.6877269 -12.6235629
##  [97] -12.5870636 -12.5742955 -12.5779435 -12.5898586 -12.6044544 -12.6217762
## [103] -12.6489228 -12.6990700 -12.7883489 -12.9315227 -13.1373601 -13.4041478
## [109] -13.7156159 -14.0383220 -14.3232947 -14.5159238 -14.5746988 -14.4899695
## [115] -14.2884211 -14.0193071 -13.7340281 -13.4719928 -13.2561370 -13.0948841
## [121] -12.9868227 -12.9259627 -12.9064607
# Check performance
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 7, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(18,1) With Weekly/Monthly Trend', AIC_val = e$aic) #ASE .000178

preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 90, lastn = T, limits = F)
eval_model(us$positivity,preds$f, model_name = 'ARMA(18,1) With Weekly/Monthly Trend', AIC_val = e$aic) #ASE .000707

rolling <- rolling_ASE(us, e, s=7, horizon=7, model_name = 'ARMA(18,1) with Weekly Trend') #ASE .000093
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000160550780936603"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.0000856798636025404"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000230595258989943"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000230459645617526"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000032800419583193"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.0000252534592547995"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## Warning: Removed 89 row(s) containing missing values (geom_path).
Signal Plus Noise

The signal plus noise model performs very well in the 7 day window, but very poorly over the 90 day window. This is a deterministic signal with a stationary mean model, and our results aren’t surprising since it doesn’t account for the weekly trend well.

preds <- fore.sigplusnoise.wge(us$positivity, max.p = 12, n.ahead = 7, limits=F)
## 
## Coefficients of Original polynomial:  
## 0.4763 0.1833 0.1282 0.1498 0.0814 -0.0857 0.0905 0.1286 -0.0125 0.0792 -0.1415 -0.1573 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9128B+0.9191B^2    1.0406+-0.0722i      0.9587       0.0110
## 1-1.3370B+0.8117B^2    0.8236+-0.7441i      0.9009       0.1169
## 1-0.3817B+0.7741B^2    0.2465+-1.1096i      0.8798       0.2152
## 1+1.3262B+0.7343B^2   -0.9031+-0.7392i      0.8569       0.3908
## 1+0.4065B+0.7162B^2   -0.2838+-1.1470i      0.8463       0.2886
## 1+1.4225B+0.5177B^2   -1.3737+-0.2106i      0.7195       0.4758
##   
## 

eval_model(us$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise',AIC_val = 0) #ASE = .000169

preds <- fore.sigplusnoise.wge(us$positivity, max.p = 12, n.ahead = 90, limits=F)
## 
## Coefficients of Original polynomial:  
## 0.4899 0.1825 0.1263 0.1449 0.0783 -0.0882 0.0842 0.1163 -0.0191 0.0749 -0.1362 -0.1288 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9009B+0.9072B^2    1.0476+-0.0687i      0.9525       0.0104
## 1-1.3188B+0.7894B^2    0.8353+-0.7543i      0.8885       0.1169
## 1-0.3724B+0.7559B^2    0.2463+-1.1235i      0.8694       0.2157
## 1+1.3275B+0.7233B^2   -0.9177+-0.7351i      0.8505       0.3925
## 1+0.4049B+0.6962B^2   -0.2908+-1.1627i      0.8344       0.2890
## 1+1.3697B+0.4724B^2   -1.4497+-0.1229i      0.6873       0.4865
##   
## 

eval_model(us$positivity,preds$f, preds$ul, preds$ll,'SigPlusNoise',AIC_val = 0) #ASE = .002678

rolling <- rolling_ASE(us, e, s=7, horizon=7, model_name = 'SigPlusNoise', model_type = 'SigPlusNoise', p = 15) #ASE .00012
## [1] "test window: 240:246, train window: 210:239"
## 
## Coefficients of Original polynomial:  
## 0.2342 -0.4855 -0.1223 -0.2171 -0.0969 -0.1119 0.3218 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.1778B+0.9139B^2    0.6444+-0.8240i      0.9560       0.1444
## 1+0.2335B+0.8023B^2   -0.1455+-1.1069i      0.8957       0.2708
## 1+1.3852B+0.6500B^2   -1.0656+-0.6349i      0.8062       0.4145
## 1-0.6751B              1.4812               0.6751       0.0000
##   
## 

## [1] "Window ASE: 0.000208585888080764"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"

## [1] "Window ASE: 0.0000342786896513529"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## 
## Coefficients of Original polynomial:  
## -0.1222 -0.3150 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+0.1222B+0.3150B^2   -0.1939+-1.7712i      0.5612       0.2674
##   
## 

## [1] "Window ASE: 0.0000312625829956592"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"

## [1] "Window ASE: 0.000268642564641623"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"

## [1] "Window ASE: 0.000057908435916914"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## 
## Coefficients of Original polynomial:  
## -0.1703 -0.2635 -0.1903 0.1825 0.0224 0.1150 0.1611 0.0453 0.0858 -0.2458 -0.4970 -0.0954 -0.3683 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.8801B+0.9434B^2    0.9965+-0.2590i      0.9713       0.0405
## 1-1.2359B+0.9403B^2    0.6572+-0.7947i      0.9697       0.1400
## 1+0.9641B             -1.0373               0.9641       0.5000
## 1+1.6980B+0.9042B^2   -0.9389+-0.4737i      0.9509       0.4256
## 1+0.9115B+0.8715B^2   -0.5229+-0.9349i      0.9336       0.3312
## 1-0.0174B+0.7523B^2    0.0116+-1.1529i      0.8674       0.2484
## 1-0.2699B+0.7264B^2    0.1857+-1.1585i      0.8523       0.2247
##   
## 

## [1] "Window ASE: 0.000120756392201658"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

Neural Network

The neural network model performs fairly well in the short term, and has the best rolling window ASE out of the bunch using minimal parameter tuning (setting lags = 7, and doing 5 fold cross validation). It is, however, the most difficult to explain, using Multilayer Perceptrons (MLP) for prediction. The example below uses two hidden nodes, both of which are autoregression lags, which apply a weigh parameter to the training input, and to each other, resulting in the prediction outcome.

ts_us <- ts(us$positivity[1:239], start = '1')
x = mlp(ts_us, lags = 7, hd.auto.type = 'cv')
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
plot(x)

preds <- predict(x, 7)
plot(preds)

eval_model(us$positivity,preds$mean, model_name = 'NNFOR', AIC_val =  0) #ASE = .00021
## Warning: attributes are not identical across measure variables;
## they will be dropped

ts_us <- ts(us$positivity[1:156], start = '1')
x = mlp(ts_us, hd.auto.type = 'cv', sel.lag = T)
plot(x)

preds <- predict(x, 90)
plot(preds)

eval_model(us$positivity,preds$mean,model_name = 'NNFOR', AIC_val =  0) #ASE = .002771
## Warning: attributes are not identical across measure variables;
## they will be dropped

rolling <- rolling_ASE(us, e, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .00009
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.0000362940126580657"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000153520775693192"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000429492839377769"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.000239811808026264"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000669272791843012"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.000169772044069679"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

Multivariate Analysis

Taking what we learned from our univariate analysis, let’s see if we can add some explanatory variables to explain more difference in the data. For Louisiana, I have selected average temperature, wind speed, precipitation, as well as significant events that I think would result in increased gatherings. These include holidays, holiday weekends, the two weeks of school start dates across Louisiana, protests resulting from the Black Lives Matter movement, and various voting days.

For the US dataset, I will only be using the significant dates data, as the weather data only makes sense when applied to a specific region, such as Louisiana.

Louisiana

MLP

la_XDF <- data.frame(
  temp = ts(xdf$AvgTempF),
  humidity = ts(xdf$AvgHumidityPercent),
  wind = ts(xdf$AvgWindSpeedMPH),
  precip <- ts(xdf$PrecipitationIN),
  gathering <- ts(xdf$GatheringBinary)
)

## 7 day forecast
ts_la <- ts(la$positivity[1:239], start = '1')
x = mlp(ts_la, hd.auto.type = 'cv', lags=7, xreg = la_XDF)
plot(x)

preds <- forecast(x, h=7, xreg=la_XDF)
plot(preds)

eval_model(la$positivity, preds$mean, model_name = 'NNFOR', AIC_val = 0) #ASE = .000664
## Warning: attributes are not identical across measure variables;
## they will be dropped

## 90 day forecast
ts_la <- ts(la$positivity[1:156], start = '1')
x = mlp(ts_la, hd.auto.type = 'cv', lags=7, xreg = la_XDF)
plot(x)

preds <- forecast(x, h=90, xreg=la_XDF)
plot(preds)

eval_model(la$positivity, preds$mean, model_name = 'NNFOR', AIC_val = 0) #ASE = .000176
## Warning: attributes are not identical across measure variables;
## they will be dropped

rolling <- rolling_ASE(la, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000535
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000971578708755374"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000125575980051967"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.000156638653876302"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.00015584632792789"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.00150668585334667"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000545042330189"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

VAR

VARselect(la$positivity[1:239], lag.max = 7, type = "both", season = 7, exogen = la_XDF[1:239,]) #AIC = -8.5337900068
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                    1             2             3             4             5
## AIC(n) -8.5264373480 -8.5337900068 -8.5269742401 -8.5200371542 -8.5124988698
## HQ(n)  -8.4425560122 -8.4439171471 -8.4311098564 -8.4181812465 -8.4046514381
## SC(n)  -8.3184445755 -8.3109406078 -8.2892682145 -8.2674745019 -8.2450795910
## FPE(n)  0.0001981888  0.0001967436  0.0001980968  0.0001994845  0.0002010039
##                    6             7
## AIC(n) -8.5091143397 -8.5027303904
## HQ(n)  -8.3952753840 -8.3828999107
## SC(n)  -8.2268384342 -8.2055978583
## FPE(n)  0.0002016965  0.0002030007
#VARselect picks p=2 using AIC and p=1 using BIC

vfit=VAR(cbind(Positivity = la$positivity, la_XDF)[1:239,], p=1,type='both', season = 7)
preds=predict(vfit,n.ahead=7)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3], 
           pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .000650

vfit=VAR(cbind(Positivity = la$positivity, la_XDF)[1:156,], p=1,type='both', season = 7)
preds=predict(vfit,n.ahead=90)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3], 
           pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .000625

rolling <- rolling_ASE(la, s=7, horizon=7, model_name = 'VAR', model_type = 'VAR', p = 1, df_XDF = la_XDF) #ASE .000281
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000469447094936274"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.000150722025793936"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000558290163363635"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000134574853250694"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.000486713927394731"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.000389727851279653"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

United States

MLP

us_XDF <- data.frame(
  gathering <- ts(xdf$GatheringBinary)
)

ts_us <- ts(us$positivity[1:239], start = '1')
x = mlp(ts_us, lags = 7, hd.auto.type = 'cv', xreg = us_XDF)
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
plot(x)

preds <- predict(x, 7)
plot(preds)

eval_model(us$positivity,preds$mean,model_name='NNFOR', AIC_val =  0) #ASE = .000204
## Warning: attributes are not identical across measure variables;
## they will be dropped

ts_us <- ts(us$positivity[1:156], start = '1')
x = mlp(ts_us, hd.auto.type = 'cv', lags =  7, xreg = us_XDF)
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).
plot(x)

preds <- predict(x, 90)
plot(preds)

eval_model(us$positivity,preds$mean, model_name='NNFOR', AIC_val =  0) #ASE = .004178
## Warning: attributes are not identical across measure variables;
## they will be dropped

rolling <- rolling_ASE(us, s=0, horizon=7, model_name = 'NNFOR', model_type = 'NNFOR') #ASE .000098
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.0000365472078758298"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 210:216, train window: 180:209"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000153293647402558"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 180:186, train window: 150:179"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000439697478267936"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 150:156, train window: 120:149"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.000256916030066162"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 120:126, train window: 90:119"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.0000639204210989819"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length
## [1] "test window: 90:96, train window: 60:89"
## Warning in preprocess(y, m, lags, keep, difforder, sel.lag, allow.det.season, :
## No inputs left in the network after pre-selection, forcing AR(1).

## [1] "Window ASE: 0.000140189117737396"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

VAR

VARselect(us$positivity[1:239], lag.max = 7, type = "both", season = 7, exogen = us_XDF[1:239,]) #AIC = -8.5191652146
## Warning in VARselect(us$positivity[1:239], lag.max = 7, type = "both", season = 7, : No column names supplied in exogen, using: exo1 , instead.
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      6      4      4      6 
## 
## $criteria
##                    1             2             3             4             5
## AIC(n) -8.3318103491 -8.4588431215 -8.4875954850 -8.5030175995 -8.4956045603
## HQ(n)  -8.2718951092 -8.3929363577 -8.4156971972 -8.4251277877 -8.4117232245
## SC(n)  -8.1832440830 -8.2954202288 -8.3093159658 -8.3098814537 -8.2876117879
## FPE(n)  0.0002407487  0.0002120323  0.0002060271  0.0002028791  0.0002043947
##                    6             7
## AIC(n) -8.5035667501 -8.4950952290
## HQ(n)  -8.4136938903 -8.3992308453
## SC(n)  -8.2807173511 -8.2573892034
## FPE(n)  0.0002027806  0.0002045136
#VARselect picks p=4 using BIC and p=6 using AIC... p=6 returns NA values, so using p=4

vfit=VAR(cbind(Positivity = us$positivity, la_XDF)[1:239,], p=4,type='both', season = 7)
preds=predict(vfit,n.ahead=7)
eval_model(la$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3], 
           pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .001218

vfit=VAR(cbind(Positivity = us$positivity, la_XDF)[1:156,], p=4,type='both', season = 7)
preds=predict(vfit,n.ahead=90)
eval_model(us$positivity, preds$fcst$Positivity[,1], pred_ul = preds$fcst$Positivity[,3], 
           pred_ll = preds$fcst$Positivity[,2], model_name = 'VAR', AIC_val = -8.5337900068) #ASE = .001590

rolling <- rolling_ASE(us, s=7, horizon=7, model_name = 'VAR', model_type = 'VAR', p = 4, df_XDF = us_XDF) #ASE .000122
## [1] "test window: 240:246, train window: 210:239"
## [1] "Window ASE: 0.000170145130090431"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 210:216, train window: 180:209"
## [1] "Window ASE: 0.0000490511908879048"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 180:186, train window: 150:179"
## [1] "Window ASE: 0.0000488062219779536"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 150:156, train window: 120:149"
## [1] "Window ASE: 0.000299162276131377"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 120:126, train window: 90:119"
## [1] "Window ASE: 0.0000907668571318613"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## [1] "test window: 90:96, train window: 60:89"
## [1] "Window ASE: 0.0000714752472492872"
## Warning in ASE$plots[x] <- plot(eval_model(data_window, pred_object, model_name
## = model_name, : number of items to replace is not a multiple of replacement
## length

## Warning: Removed 89 row(s) containing missing values (geom_path).

Model Comparison

Type Model Rolling ASE Dataset
Univariate Louisiana AR(15) Weekly .000277
Univariate Louisiana ARIMA(14,0,1) Weekly .000679
Univariate Louisiana SigPlusNoise .000418
Univariate Louisiana MLP .000637
Multivariate Louisiana MLP .000601
Multivariate Louisiana VAR .000281
Univariate National AR(9,1) .000096
Univariate National ARIMA(18,0,1) Weekly .000093
Univariate National SigPlusNoise .000120
Univariate National MLP .000090
Multivariate National MLP .000098
Multivariate National VAR .000122

Final Model Selection

LA model: Weekly AR(15)

𝜑(1−.682𝐵−.23〖3𝐵〗2−.179𝐵3+.144𝐵4+.035𝐵5−.073𝐵6+.644𝐵7−.364𝐵8−.082𝐵9+.066𝐵10+.128𝐵11−.015𝐵12−.216𝐵13+.341𝐵14−.306𝐵15)(1−𝐵7)=𝑎𝑡 - 7-day ASE = . 000647 - 90-day ASE = . 000202 - Rolling ASE = .000277 - AIC = -8.488985

US model: AR(9,1)

𝜑(1−1.271𝐵+.163𝐵2+.070𝐵3−.066𝐵4−.003𝐵5+.164𝐵6−.161𝐵7−.114𝐵8+.229𝐵9)=(1−.882𝐵)𝑎𝑡 - 7-day ASE = .000121 - 90-day ASE = .001150 - Rolling ASE = .000096 - AIC = –8.625600

Summary

In both cases, the Univariate models outperformed their multivariate counterparts. This may suggest that I selected poorly correlated features to include in the model. But as a testament to the time series algorithms, they perform just as well as the much more complicated neural networks.

e <- est.arma.wge(la.s7, p=15, q=0)
## 
## Coefficients of Original polynomial:  
## 0.6819 0.2329 0.1795 -0.1438 0.0353 0.0725 -0.6438 0.3639 0.0825 -0.0664 -0.1285 0.0148 0.2163 -0.3415 0.3063 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.8995B+0.9868B^2    0.9625+-0.2950i      0.9934       0.0473
## 1+0.9661B+0.9242B^2   -0.5226+-0.8994i      0.9614       0.3338
## 1+1.3916B+0.9115B^2   -0.7634+-0.7172i      0.9547       0.3800
## 1+1.8666B+0.8970B^2   -1.0405+-0.1795i      0.9471       0.4728
## 1-0.9455B              1.0576               0.9455       0.0000
## 1-0.1935B+0.8186B^2    0.1182+-1.0989i      0.9048       0.2329
## 1-1.3104B+0.7391B^2    0.8865+-0.7531i      0.8597       0.1121
## 1-0.5571B+0.7181B^2    0.3879+-1.1145i      0.8474       0.1967
##   
## 
preds <- fore.aruma.wge(la$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 30, lastn = F, limits = F)

ggplot() + 
  geom_line(aes(us$date, us$positivity), color='#F8766D') +
  geom_line(aes(
    seq(as.Date("2020/11/12"), by = "day", length.out = 30), 
    preds$f), color='#00BFC4') + labs(title = 'Lousiana 30 day predictions', y = 'Positivity Rate', x='Date')

e <- est.arma.wge(us$positivity, p=9, q=1)
## 
## Coefficients of Original polynomial:  
## 1.2711 -0.1631 -0.0699 0.0662 0.0028 -0.1636 0.1612 0.1144 -0.2295 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.9861B+0.9886B^2    1.0045+-0.0501i      0.9943       0.0079
## 1-0.0326B+0.7214B^2    0.0226+-1.1772i      0.8493       0.2469
## 1-1.1925B+0.6943B^2    0.8588+-0.8384i      0.8332       0.1231
## 1+1.2068B+0.6319B^2   -0.9548+-0.8190i      0.7949       0.3872
## 1+0.7333B             -1.3636               0.7333       0.5000
##   
## 
preds <- fore.aruma.wge(us$positivity, phi = e$phi, theta = e$theta, s=7, n.ahead = 30, lastn = F, limits = F)

ggplot() + 
  geom_line(aes(us$date, us$positivity), color='#F8766D') +
  geom_line(aes(
    seq(as.Date("2020/11/12"), by = "day", length.out = 30), 
    preds$f), color='#00BFC4')+ labs(title = 'US 30 day predictions', y = 'Positivity Rate', x='Date')